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Introduction 


This unit is devoted to the study of numerical methods for finding approximate 
solutions of differential equations. Only first-order differential equations of the 
form 


y= m(x,y), 


with (xo) given, will be considered but the methods in this unit can be adapted to 
find approximate solutions of higher-order differential equations. (A differential 
equation of any order may be treated as a system of first-order differential 
equations, as described in Section 5 of Unit 6.) 


The primary objective of this unit is to enable you to use numerical methods 
intelligently to obtain approximate solutions. To do this it is necessary to study 
several different methods 

(i) to see how they work, 

(ii) to determine whether one method is better than another, 

(iii) to see when things go wrong. 


All the methods in this text lead to recurrence relations which can be used step by 
step to obtain an approximate solution at equally spaced points Xo, X1, X2,+.. « 
As in Unit 18 we denote the values of this approximate solution by Yo, Yi, Yo... « 
It is important to compare this solution with the true solution y = y(x). In 
particular we are interested in the difference between the approximate value, ¥,, 
and the true value, y,, at x). 


Study guide 


This unit is intended to be read sequentially. In Sections | and 2 we derive some 
numerical methods and use them to solve simple problems. The television section 
(Section 2) illustrates some of these methods and how they can be used to solve 
the logistic equation, Before viewing the programme it is advisable to have read 
Subsection 1.1 on the trapezoidal method and also the pre-broadcast notes. 


In Section 3 and in the tape session in Section 4 we look at the methods of the 
first two sections to discover why they work, when they go wrong and what sort of 
accuracy you could expect. 


In the remainder of Section 4 we look at Simpson's method which is more 
accurate than most of the one-step methods but can be difficult to use. 


Section 5 describes a computer package which you can use to solve differential 
equations numerically. One of the TMA questions for this unit may well require 
you to use this package. 


1 Numerical methods 
1.1 Integration methods 


In Unit 18 you met three integration methods for evaluating the integral 
rt 
l= S(x)dx. 


xp 


They were Euler's method, the trapezoidal method and Simpson's method. In this 
section we will see how the first two of these methods can be used to solve first- 
order differential equations of the form 


y' = m(x,y). (1) 
In Section 4 we will examine Simpson's method. 


If we integrate both sides of Equation (1) between x, and x,,, we have 


[vax < [mena @) 


You can use the computer 
package at Summer School. 


‘4 dy 
y’ is short for = 
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The left-hand side of this equation can be integrated to give 

fry dX = Ven —Ye y, is short for y(x,). 
Hence we ny write Equation (2) as 

Vert Yr = fF messy (G) 


Thus, if we can determine the integral on the right-hand side and if we know Vrs 
we can compute y,,, using Equation (3). The difficulty with this approach is that 
we do not know the function y, so it is unlikely that we will be able to determine 
the integral. 


For example, the differential equation 
y =xy with y(0) =0 

has m(x,y) = xy. 

Equation (3) becomes 


Xe 


1 
xydx, 


Since we do not know y(x) we cannot proceed to determine an exact recurrence 

relation, but we can determine an approximate one if we approximate the integral. 

In the last unit we saw that the simplest method for approximating integrals was 

Euler’s method (see Figure 1) giving fi 


[ferae= nen 


where h = x,4, — x,. Kix) 


In the same way we can approximate the integral in Equation (3) using Euler’s 
method as 


3, are 


rt b—h—4 
f m(x,y)dx = hm(x,,),). 

Ne Figure 1. The shaded area is 
With this approximation Equation (3) becomes ian Ly eee rue 


Yea. = ¥, + hm(x,, ¥,). 


(I have used capital letters, Y, and Y,,,, here because Y, will only be an 
approximation to the true solution.) But this recurrence relation is just Euler's 
method for solving a first-order differential equation, as described in Units 2 and 
18. This new derivation illustrates how integration methods can be used to solve 
differential equations. 


To simplify the notation, we define 
Y, = m(x,, ¥,) 

and write the recurrence relation for Euler's method as 
Y+1=¥, + AY}. 
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Euler’s method 
1. You are given a differential equation 
y= m(x,y) 
and a value for y(xo). 


2. To apply Euler’s method, choose a step-size, h, and calculate 
Y,, Y... from 


Yi =¥,+hY, 


where Y{ = m(x,, Y,), 
Yo = y(Xo), 
and X,=Xo trh. 


3. Y, is an approximation to y,. 


Another integration method is the trapezoidal method. From the results in Unit 18 
it is clear that this method for approximating the integral by a trapezium (see 
Figure 2) as 


eet h ; 
[Fea = F004) + For) 
Xr 
is generally more accurate than Euler’s method. We would thus expect to get more 
accurate results if we use this method to solve differential equations. 


If we apply this method to approximate the integral in Equation (3) we have 
[7 missonas = Himes.) + M(Xp415Vr+ 1) 
This gives rise to the recurrence relation 
Yon = Ye + ElmG.p, Y) + mbsreas Fra] 
or, more simply, 
Yer = Yet Ave + b Grn (4) 


However there is a difficulty in using the trapezoidal method because Equation (4) 
contains a term in Y4, (= m(x,+1, Y,+1)) on the right-hand side, We do not 
know this value until we have computed Y,,,. Equation (4) would have to be 
solved for Y,,,. For linear equations this can be done fairly easily, as Example 1 
shows. 


Example 1 
Use the trapezoidal method with h = 0.1 to determine the approximate solution 
forx,=rh, O0<r<10, of 


with y(0) = 0. 


Compare the answers with the true solution, y = e* — x — 1. 


yoxty 


Solution 
We want to apply the recurrence relation (4). The differential equation tells us that 


m(x,y) =x + y 


andso Yj; x, + Y, 


and Yie1 =Xr41 + Vers 


We substitute for Y and Y;,, in Equation (4) to get 


h 
Your = Ye + 5 [ler + Ye) + Oras + Yeo): 


37 Xr 


Fh 


Figure 2. The shaded area is 
the trapezoidal approximation 
to the area under the curve. 


7% 
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To solve for Y,,, we collect up the terms in Y,,, and Y, as 
h h\ oh 
Y(t -4) = v(t +4) +5 (%e + Xr41): 


. _ (CL + h/2)¥, + (h/2)x- + Xr41) 
Le. Yor= 1-h2 


We can use this recurrence relation to generate the approximate solution. We have 
(i) h=0., 

(ii) x,=X9 +rh=O0.1r, 

We can simplify the recurrence relation to 


1,05Y, + 0.05(0.1r + 0.1(r + 1)) 


Yer = 095 


7a 2 te +1 


The computation is now straightforward, starting with Y) = 0. The following table 
shows how the approximate solution, generated using the recurrence relation for 
the trapezoidal method, compares with the true solution. The results for Euler's 
method and the Taylor series method of order 2 are also given for comparison. 
(These were obtained in Unit 18.) 


Solutions to y’ = x + y with y(0) = 0, h = 0.1, 


trapezoidal true Taylor series method 

| ee method solution | Euler’s method of order 2 
0 | 0.0 0.0000 0,0000 0,0000 0.0000 
1/01 0.0053 0.0052 0.0000 0.0050 
2) 02 0.0216 0.0214 0.0100 0.0210 
3 | 03 0.0502 0.0499 0.0310 0.0492 
4) 04 0.0923 0.0918 0.0641 0.0909 

5] 05 0.1494 0.1487 0.1105 0.1474 

6 | 0.6 0.2230 0.2221 0.1716 0.2204 

7) 07 0.3149 0.3138 0.2487 0.3116 

8 | 08 0.4270 0.4255 0.3436 0.4228 
9} 09 0.5615 0.5596 0.4579 0.5562 
10 | 1.0 0.7206 0.7183 0.5937 0.7141 


By comparing columns 3 and 4 we see that the results for the trapezoidal method 
are reasonably accurate although, as expected, the error increases as x increases, 
The approximate value of y at x = 1 is correct to two decimal places, Now 
compare the true solution with the results obtained using Euler’s method, For this 
problem Euler’s method gives much worse results and a much smaller step-size 
would be needed to obtain reasonable accuracy. On the other hand the results for 
the Taylor series method of order 2 are about as accurate as the results for the 
trapezoidal method. 


The manipulation required to use the trapezoidal method to solve the general 
linear differential equation of the form 


y’ = Ixy + K(x), (5) 
using a step-size h, can be deduced in the same way as in Example 1. The 
recurrence relation (4) for the trapezoidal method is 

h 

Yas = Ye +5(¥e + Yea). 6) 
Using Equation (5) we have 

Y=1,Y, +k, 
and Yes = Trea Yeva + heat 


This can be expressed as 


21 iz 1 
Yrs =i9%r + 95 + I99 


which is a linear recurrence 
relation. 


1, is short for I(x,), k, is short 
for k(x,) and so on, 
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Substituting for Y; and Y/,, in Equation (6) gives 
h 
Yor =H +5 + be) + (bes Yeor + hrv)] 


which can be solved for Y,, as 


y., < LEM TYs + (h/2We + heros) 
a 1=(/2)ra1 ; 


The trapezoidal method for linear differential equations 
1. You are given a linear differential equation 
y’ =I(x)y + k(x) 
and a value for y(Xo). 


2. To apply the trapezoidal method, choose a step size h and 
calculate Y\, ¥2,... from 


_ Ut (h/2)b VY, + (h/2Mke + ra) 
By 1 = (h/2)leva 
where Yo = y(xo), and x, = Xo + rh, 


Yi+a 


3. Y, is an approximation to y,. 


We have seen two methods in this section: Euler’s method, in which the recurrence 
relation can be applied directly, and the trapezoidal method, which is difficult to 
use except for linear equations. The difference between the two methods, which 
makes the trapezoidal method harder to use, can be described by saying that the 
trapezoidal method is implicit whereas Euler's is explicit. 


In general any numerical method for solving differential equations in which the 
recurrence relation involves Y;,, or higher derivatives at x,,, is called implicit, 
This usually implies that we cannot write down the value of Y,,, directly. On the 
other hand a method is called explicit if ¥;,, or higher derivatives at x,,, do not 
appear in the recurrence relation for Y,,,. In this case Y,,, is defined explicitly in 
terms of previous values. 


Exercise 1 
Use the trapezoidal method with h = 0.2 to approximate the solution to the differential 
equation 


y' =3y + sinx with y(0) = 0, 
at x = 0,2, 0.4 and 0.6, 
The correct solution is 


which gives (0.2) = 0.02460, y(0.4) = 0.12308 and (0.6) = 0.35304. 
Compare the approximate solution with the true solution by calculating the errors. 


Exercise 2 
Write down the formula to be used when the trapezoidal method with h = 0.1 is to be 
applied to the differential equation 

y' =siny with y(0) = 1. 
What is the difficulty with this method if you try to use the recurrence relation to generate 
the approximate solution? 
[Solutions on p. 50] 


1.2 Local truncation errors 


So far in this course you have met the Taylor series methods (of which Euler’s 
method is one) and the integration methods (of which Euler’s method is again 
one). By the end of the unit you will be familiar with two more methods. With so 
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many methods available what is the best method to use? This question is by no 
means easy to answer and I shall devote most of Section 3 to examining some of 
the pitfalls. For the moment we will assume that we just have the choice between 
the Taylor series methods and the trapezoidal method and want to make some 
simple comparison between them. 


Ideally for a given problem we would like to compare the error in the nth term in 
the sequence of approximations generated by method A with the error in the same 
term given by method B. 


Yi 


fe 


This is the error 
in the computed 
solution after n steps 


of: 


°% The point (x), y,) has, for 
simplicity, been labelled y, 


and so on. 


° 
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Figure 3. Graph of the approximate solution and the true solution. 


However, to do this we would need to know the true solution, which is not 
normally available, so we need an alternative method of comparison. There is a 
fairly simple yardstick by which we can judge a numerical method and that is to 
estimate the error introduced in one single step; if we use a method which 
introduces a comparatively small error at each step then the overall error is likely 
to be comparatively small. 


From Unit 2 Section 1 on direction fields we know that there is a whole family of 
solution curves. Thus after r steps we can assume that our approximation Y, lies 
On some solution curve (not necessarily the one defined by the initial condition), 
The error we are going to measure is the distance of Y,,, from this particular 
solution curve (see Figure 4). 


To illustrate the technique we look at Euler’s method, given by 
Yo1=Y¥,+ hyp. 


Let y = (x) be the solution curve in Figure 4 on which Y, lies. We have 


® Y=», &, Her 
(i) ¥p = m(x,, ¥,) = M(x 94) = Yee Figure 4, Graph showing the 
Thus error in a single step. 
Yer =yrt hy 7”) 
The true value, y,.,, can be found by using the Taylor series expansion at x, to 
give 
h? he 
Verve Te te Meck ve Vr ht 5 (8) 
2 6 
The error in Y,,, is given by ¥,,, —y,4, which we can write, using Equations (7) 
and (8), as 
h? BB 
Teer Yeo = Oe + hye) — [ye + hye + Tye + Eye to 
he 3 
ee 9 
Qo" ey (9) 


This error in Y,,, is called the local truncation error because it is 


(i) local in the sense that we are looking at the error in a single step, i.e. in the 
locality of x,; 
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(ii) a truncation error in the following sense: by comparing Equations (7) and (8) 
it can be seen that we could have obtained the value of Y,,, by truncating the 
Taylor series in Equation (8), 


Look at the local truncation error for Euler's method (Equation (9)). The first 
term in this series, (—h?/2)y’, is the most important because it contains the lowest 
power of h. For sufficiently small values of h all the other terms become relatively 
small compared with this first term. For this reason we give the first term a special 
name and call it the principal term in the local truncation error. 


When comparing two methods it is the principal terms in their local truncation 
errors which will normally be compared and in particular the power of h in the 
principal term. The method with the higher power of h in its principal term will 
give a smaller local truncation error if h is sufficiently small. 


I must emphasize that the local truncation error does not, in itself, give us any 
idea of the accumulated error after n steps. We shall deduce some results for 
accumulated errors in Section 3, The local truncation error does give us some 
practical information concerning the choice of step-size. Suppose we halve the 
step-size, h, in Euler’s method. Since the principal term involves h?, halving the 
step-size reduces the local truncation error by a factor of 4, approximately, 
Remember, however, that halving the step-size means that we have to do twice as 
many computations to determine the approximation at a given value of x, so that 
there is more accumulation of errors. 


In the following example we shall see how the accumulated error behaves for two 
different values of h. 


Example 2 
Apply Euler's method to the differential equation 
y=xty with y(0) = 0, 


using step-sizes h = 0.1 and h = 0.05 to compute the values of Y, at x = 0.1, 0.2 
and 0.3. 


Solution 
h=01 h=0.05 
| r Xr i; r Xy. | ¥, 
0 ) 0 0 0 0 
1 0.05 0.0000 
1 0.1 0.00 2 0.1 0.0025 
3 0.15 0.00763 
2 0.2 0.01 4 0.2 0.01551 
5 0.25 0.02628 
3 0.3 0.031 6 0.3 0.04010 


The true solution at x = 0,3 is y = 0.0499. The graph in Figure 5 opposite 
indicates the growth of errors, 


Exercise 3 
The Taylor series method of order 2 is given by 


ree 
Yay = Y+hY, + >Y¥r. 


Calculate the local truncation error for this method. Write down the principal term in the 
local truncation error. 
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0 On 0.2 0.3 x 


Figure 5. Graph comparing the numerical results for two 
values of h in Example 2, 


Exercise 4 

Calculate the local truncation error for the Taylor series method of order 3. Write down the 
principal term in this error, By what factor would you expect the local truncation error to 
decrease if the step-size were halved? 


Exercise 5 
Calculate the local truncation error for the Taylor series method of order n, 
(The result of this exercise is important and is repeated in the summary.) 


Exercise 6 
The Taylor series method of order 2 was applied to the differential equation 
yloxt+y with y(0) = 0. 


Two sets of results were obtained for step sizes h = 0.1 and h = 0.05 for values of x up to 
x=03. 


h=01 h= 0.05 
| r x y, r x % 
0 0 0 0 0 0 

1 0.05 0.00125 

1 0.1 0.005 2 O41 0.005127 

| O15 0.011764 

2 0.2 0.021025 4 0.2 0.021305 

5 0.25 0.033897 

3 03 0.049232 6 03 0.049696 


‘The true values are y = 0.005171 at x 
y = 0.021403 at x 
y = 0.049859 at x = 


Compute the errors at x = 0.1, 0.2 and 0.3 for the two step-sizes. By approximately what 
factor have the errors been reduced by halving the step-size? 


[Solutions on p. 51] 


Although halving the step-size 
reduces the error per step by a 
quarter (approximately), the 
error in the value of ¥ at 

x = 0,3 has only been 

reduced by 

(Cees tely), We 

shall investigate this further in 
Section 3. 
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1.3 Local truncation errors for implicit methods 


Local truncation errors for implicit methods, such as the trapezoidal method, 
cannot be calculated using the techniques you have met in the last subsection. In 
the trapezoidal method we have 


Y= 


AE + Yee (10) 
If we assume that Y, lies on a solution curve we have 
h 
Yor = Yet 50e + Yr+1) (11) 


and the difficulty is that we do not know enough about Y;,, to be able to 
compare this formula for ¥,,, with the Taylor series expansion for y,,,. However 
we can compute the principal term in the local truncation error by replacing ¥;,, 
by y/41 in Equation (11), Although the error in Y;,, affects the other terms in the 
local truncation error, it does not affect the first term since Y},, is multiplied by h 
in Equation (11), In effect we are saying that the principal term in the local 
truncation error can be computed by assuming that, in Equation (10), all the 
right-hand side values can be replaced by their true values. 


The Taylor series expansion for y;,,, is 


h? 
Voor = Yet hye + oye toe 


Substituting this in the approximate expression for Y,,, gives 


h 
Yor = Yet 50 + roi) 


AOR atnaptay costae 
Hot Bet Ue + hye + aye + °°) (12) 
Now the true value for y,,, is given by the Taylor series at x, as 
ae 
Yess =Vet Wr + Vet Eye toe (13) 


Subtracting Equation (13) from Equation (12) gives 


hen renee ee fis nioerge ue 
Yous — Yorn = Ve tae + ahve + hye + Ve} — Ye — hye — SVE — Ee — 


hs 
mae ta 


Thus the principal term in the local truncation error for the trapezoidal method is 
(h?/12)yp. 


The following procedure summarizes the technique for computing the principal 
term in the local truncation error and can be used for both implicit and explicit 
methods. For an explicit method you get all the terms in the local truncation error 
whereas for an implicit method you can only get the principal term. 
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To find the principal term in the local truncation error 


1. Write down the recurrence relation for the method as a 
formula for ¥,,,. 


2. Replace Y by y in all symbols on the right-hand side of this 
recurrence relation. 


3. Expand this new recurrence relation as a power series in h 
using Taylor series about x,. This gives Y,,, in terms of y,, 
Verve 

4. By subtracting the Taylor series 

2 


agllana 
Vea = Vet hye tye te 


express ¥,,,; — y,4, in powers of h, 


5. The term in the lowest power of h is the principal term in the 
local truncation error. 


We can compare, for example, the trapezoidal method with Euler’s method. Since 
the principal terms in the local truncation error involve h? and h? respectively we 
would expect much better results from the trapezoidal method for sufficiently 
small values of h. 


Exercise 7 
Consider a method which uses the recurrence relation 
Yor = Ye thYeos. 


Determine the principal term in the local truncation error for this implicit method. What 
other method have you seen whose principal term involves the same power of h? 


Exercise 8 


(i) For the trapezoidal method, if you halve the step-size what is the approximate factor 
by which the error per step reduces? Which Taylor series method would reduce the 
error per step by the same factor? 


(ii) Using the results of Exercise 6, by what factor would you expect the accumulated errors 
to reduce if the step-size were halved for the trapezoidal method? 


Exercise 9 
The trapezoidal method was applied to the differential equation 
yax+y with y(0) = 0, 


Two step-sizes, h = 0.1 and h = 0.05, were used as in Exercise 6. The following tables of 
results were obtained. 


h=01 h = 0.05 
r SS y, r x y, 
eit 
0 0 i) 0 ) 0 
1 0.05 0.001282 
1 4 0.005263 2 Ou 0,005194 
3 0.15 0.011871 
2 02 0.021607 4 02 0.021454 
5 0.25 0.034092 
3 03 0.050197 6 03 0.049943 


Comparing these results with the true solution at x = 0.1, 0.2 and 0.3 we have the following 
table of errors. 


xX, h=01 h=0.05 
01 0.000092 0.000023 
0.2 0.000204 0.000051 


03 0.000338 0.000084 
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By what factor (approximately) have the accumulated errors been reduced by halving the 
step-size? 
[Solutions on pp. 51-52] 


Summary of Section 1 


1. Integration methods for the numerical solution of differential equations are 
derived from numerical integration techniques, 


2. The trapezoidal method is an integration method and gives rise to the 
recurrence relation 


h 
Yo. = Yt Vet Yea) 


3. The trapezoidal method can be applied to the linear differential equation 
y =Ix)y + k(x) 
using the procedure on page 8. 


4. The trapezoidal method is an example of an implicit method. An implicit 
method for solving a differential equation gives rise to a recurrence relation 
involving Y},, or higher derivatives at x,,,. A method is explicit if the 
recurrence relation does not involve Y;,, or higher derivatives at x,,,. Implicit 
methods are difficult to use with non-linear equations. 


5, The local truncation error for a numerical method for solving differential 
equations is the error induced in Y,,, assuming that Y, lies on some solution 
curve, y, so that Y, = y,. 


6. The principal term in the local truncation error is the term with the lowest power 
of h. The principal term can be determined using the procedure on page 13, 


7. The principal terms in the local truncation error for various methods are as 
follows: 


method principal term 
h? 
(i) Euler’s method meyes 
3 
(ii) Taylor series method of order 2 = r aie 
- ; nee , 
(iii) Taylor series method of order n am era 
(iv) Trapezoidal method re 


2 The predictor-corrector method 
(Television Section) 


2.1 Pre-broadcast notes 


The television programme for this unit investigates ways of obtaining numerical 
solutions to the logistic equation 


, y ‘ 
= i = 100. 
y 105 ial with y(0) = 100, 
This equation models the growth in a population where the initial population is 
100 and the population stops growing if it reaches 1000. 


This equation has an analytic solution which you will be asked to find in 
Exercise 1, so you may be wondering why we should want to determine numerical 
approximations for this problem. There are several reasons: 


You met this 
Unit 3 Section 


uation in 
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(i) We can compare our numerical solution with the true solution. 


(ii) If we have to obtain numerical approximations to a closely related 
differential equation which cannot be solved analytically, we may be able to 
use information we derived solving this one. 


(iii) The non-linearity in this equation will highlight some of the difficulties in 
using numerical methods and how these difficulties can be overcome. 


The television programme is divided into three parts: 


(i) The integration methods developed in Subsection 1.1 will be reviewed and 
the methods applied to the logistic equation. An important feature of the 
programme is the way in which the graphs of y' and y inter-relate. 


(ii) The failure of the trapezoidal method gives rise to a method based on using 
Euler’s method and the trapezoidal method together in a method known as a 
predictor-corrector method, We will look at this method in more detail after 
the programme. 


(iii) The final part of the programme looks at the importance of choosing a small 
enough step-size, h, so that the numerical solution is qualitatively the same 
as the true solution. This topic will be explored more fully in Section 3. 


Exercise 1 
Use the method of separation of variables to solve the differential equation 
y 
y = 1011-2}, 
2 »( 1000) 


What is the particular solution if »(0) = 100? 


Sketch a graph of the particular solution for 0 < x < 1. (The easiest way to do this is to 
plot a few points.) 


[Solution on p. 52] 


Now watch the television programme ‘Integrating by numbers’. 


2.2 Programme notes 


Euler's method, as an integration method, can be considered as a means of 
building up a graph of the approximate solution, y, by calculating the approximate 
area under the graph of y’, as described in Subsection 1.1. In Subsection 1.1 we 
derived the basic equation for an integration method as 


Vs =Yor ij yy dx (from Equation (3) of Subsection 1.1) 


where y’ = m(x,y) is the differential equation to be solved. 


Euler's method was derived by approximating the integral by the area of a 
rectangle of width h = x; — Xo and height y% (see Figure 1(a)), 


x 
ie. f y'dx = hyo. 
Xo 
The value of yo can be found by substituting the initial condition yo in the 
differential equation to give 
Yo = m(Xo, Yo). 


The value of the area, hy} is added onto the value of yp to give an approximation 
for y; so that 


Yi = Yo + hyo 
as shown in Figure 1(b). 


Once we have found Y; we can repeat the process to find ¥{, then Y, and so on. 
In this way we can build up the graphs of both y and y’. Summarizing the method 
we have a three-stage procedure: 


Figure 1(a). Graph of y’ 
against x showing the area for 
Euler’s method, 


* 
hg 


Figure 1(b). How yo is 
incremented to obtain Y,. 
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(i) Evaluate Y; as Y; = m(x,, Y,) and plot this point on the y’ graph. 
(ii) Calculate the area of the rectangle height Y; and width h = x,,, — x,. 


(iii) Increment Y, by the value of this area to obtain Y,,,. 


yh 


Stage (iii) 
i +t 
Y, 


Figure 2. Schematic diagram for Euler's method. 


If we try to get the same geometrical insight into the trapezoidal method we come 
unstuck because the method requires both yg and Yj in order to be able to 
calculate Y, using the formula 


h 
= yo + 5(vo + Yi). 


To find Y, we must increment yo by the area of a trapezium and we cannot find 
the area of the required trapezium geometrically because the only information we 
have is the values of yo and yj = m(xo, Yo). 


For linear equations we can do some algebraic manipulation to get round this 
difficulty, as described in Subsection 1.1. For the logistic equation we could do 
some algebraic manipulation and determine Y, by solving a quadratic equation. 
This process is much more involved than Euler's method and in the programme 
this method is not pursued. 


The second part of the programme indicates how we can get round the difficulty 
in using the trapezoidal method for non-linear differential equations by using a 
combination of Euler's method and the trapezoidal method. To apply the 
trapezoidal method we need Y{ in order to compute Y, as 


h 
Y= yo +506 + Yi) 


and Y; depends on Yj. The way round this difficulty is to use Euler's method to 
obtain a crude approximation for Y{ which we can then use in the trapezoidal 
formula. Because this crude approximation is multiplied by h in the formula we 
can still get quite a good approximation, Y,. 


Example 1 
To see how the method works I will compute a value for Y, at x = 0.1 for the 
logistic equation using this new method. 


The important thing here is to 
see how the graphs of y and y' 
inter-relate. 
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Solution 


The method is in four stages: 


(i) We know that Yj = 100, so I can evaluate Yj using the a 
differential ti ¥g= 10% (1 XY oe Y= (xy. Yo) 
ifierenti Ui Ss = — meee 
equation as Y¢ Aj 1000 giving te 


Yj = 10% 100{1 — ioe = 900. | me | 

Y{ enables us to calculate the area of the rectangle in Figure i 

3(a). w=0 4 =O1 x 
Figure 3(a) 

(ii) Using Euler's method I get a crude approximation for Y, in v4 
Figure 3(b) as Grom) 
Yt = Y +hY¥o = 100+ 0.1 x 900 
= 190, ra 
I've labelled this with an asterisk to indicate that it is a A ieee 
crude approximation. We are going to do better. ———— 
0 or x 

Figure 3(b) 


(iii) This is not a particularly good approximation, However I 
can use it to obtain a value for the derivative Y’, using the —_ en ye 
1 YN) 


logistic equation, as (Yt! = man. 
ean: a 
r= 1or4(t - a - 10 x 190(1 - ionik _ 
= 1539. 
This is shown in Figure 3(c). ‘iy is 
Figure 3(c) 


(iv) Look at Figure 3(d). Since we now have values for Yj; and 
Yf’ we can calculate the area of the trapezium. Adding it to vy 
Yo we obtain an improved approximation Y, as 


watt Yi 
; 
; 
i 


%=%+ 3( + rr} 


ea 
01 } 
100 + (900 + 1539) ; 


0 On 
EEE, Figure 3(d) 
as shown in Figure 3(e). 
This is much closer than Y¥ to the true value of y at Ya 
x = 0.1, which is 231.97 correct to 2 decimal places. M=% 
Y, 
The method outlined above is an example of a predictor-corrector method. In : 
Example 1 Euler's method was used to predict a value for Y,. This value for Y, : 
was then corrected using the trapezoidal method. 7 
Predictor-corrector methods provide one of the most widely used techniques for ek ae 
solving differential equations. They are based on the idea of using an explicit & ae) 
igure 3(e 


method to obtain a prediction and then using a more accurate implicit method to 
correct this prediction. 


The four stages of the procedure I used in Example 1 can be summarized as: 
(i) Evaluate Y¢ using the differential equation. 

(ii) Predict Yf using Euler’s method. 

(iii) Evaluate Y}' using the differential equation. 


(iv) Correct Y, using the trapezoidal method. 
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Here is the same procedure in diagramatic form. 


(i) EVALUATE (iii) EVALUATE 


Area is 
AY 


xy oN xy XX 


Yi = my. YP) 


mx. Yo) 


(ii) PREDICT 
y, 


xo ‘ x x xy x 
EULER TRAPEZOIDAL 
Y¥}= Yuth¥i ¥) = Yo +4 (Yo ¥19) 


Figure 4. Schematic diagram of the Euler-trapezoidal predictor-corrector method, 


Example 2 

When the procedure above was used to compute Y, for the logistic equation, the 
value Y, = 221.95 was obtained. I can now use this value as the starting point in 
the computation of ¥,. 


Solution 
We can use the four stages again as follows: 


: hs oa 221.95 
(i) Evaluate Y; = toy (1 — 7000} = 10 x 221951 000, 
= 1726.882. 
(ii) Predict Y$ using Euler’s method as 
Ys=Y¥, +hYj hY; is the area of the 
rectangle, 

= 221.95 + O.1 x 1726.882 

= 394.6382. 
* ae vee 394.6382 
(iii) Evaluate Y}' = iors (1 7000] = 10 x 39463821 7000 


= 2388.9889. 
(iv) Correct Y, using the trapezoidal method. 
h a ny 4 
y%=% +844 yt’) 3 (Vi + YP’) is the area of the 
2 trapezium. 
= 221.95 + 0.05(1726.882 + 2388.9889) 
= 427.74355. 
Exercise 2 
Follow the method of Example 2 to compute values for Y; and Y, for the logistic equation. 
[Solution on p. 52] 


Thus we can use the Euler-trapezoidal predictor-corrector method to compute Y,, 
Yp,.--, Yq) a8 summarized in the following procedure. 
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The predictor-corrector method (Euler-trapezoidal) 


1. You are given a differential equation 


y= m(x,y) 


and a value for y(xo). 


2. To apply the predictor-corrector method, choose a step-size h, 
and calculate Y,, Y2, Y3,... by repeating the following cycle of 
steps starting with Yy = y(xo): 


(i) Evaluate Y; = m(x,, Y,). 
(ii) Predict using Euler’s method: 


Yt. = ¥, +hY}. 


(iii) Evaluate Y}4, = m(x,, Y¥+1). 


(iv) Correct using the trapezoidal method: 


h 
Yer = Year+ Yt). 


3. Y, is an approximation to y,. 
ee 


Note that any predictor-corrector method is an explicit method, (You will be 


asked to show this for the Euler-trapezoidal method in Exercise 5.) 


In general, predictor-corrector methods use an explicit method to predict Y¥,, 
and a more accurate implicit method to correct Y,,;. 


Although the derivation of the local truncation error for a predictor-corrector 
method is beyond the scope of this course, it can be shown that the principal term 
in the local truncation error for the Euler-trapezoidal method involves h® (as does 
the principal term in the local truncation error for the trapezoidal method). 


The following table gives the computed values for the logistic equation for each of 
the three methods discussed in the programme and also the true solution, obtained 
using the formula from Exercise 1, The results for the trapezoidal method were 

obtained by solving a quadratic equation at each step. 


Solutions to the logistic equation, with h = 0.1, 


Xr 


Euler 


true solution 


trapezoidal | predictor-corrector 
100 100 
235 222 
448 428 
681 660 
852 823 
942 911 
980 956 
993 978 
998 989 
999 994 
1000 997 
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T have plotted these points on a graph shown in Figure 5. Note that the 
trapezoidal method gives the best results with the predictor-corrector method 
initially giving much better results than Euler’s method. 


0 On 0.2 0.3 04 0.5 0.6 0.7 0.8 09 10 x 
Figure 5. Comparison of results for the logistic equation. 
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In the final part of the programme we look at what happens if we use different 
step-sizes, h. The first method we look at is Euler’s method. The four graphs in 
Figure 6 show the results for h = 0.1, 0.15, 0.2 and 0.25. The logistic curve is also 
shown for comparison. 


th 


0 04 08 1.2 16 20 x 0 05 10 Ls 2.0 25 x 


Figure 6. Euler's method for different step-sizes, 


The results show that, for h = 0.1 and 0.15, the approximate solution at least has 
approximately the same qualitative behaviour as the analytical solution, The saw- 
tooth behaviour near y = 1000 of the solution for h = 0.2 and 0.25 is clearly 
undesirable. This will be discussed in Section 3 but the graphs do indicate that, for 
some problems at least, we have to be careful about our choice of step-size if we 
want to get reasonable results. 
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The trapezoidal method, on the other hand, behaves well for this problem as the 
graphs in Figure 7, with step-sizes h = 0.1, 0.2, 0.5 and 1.0, illustrate. 


FR 


a8: 


Seii eseenaael atsttaees 


0 0.5 1.0 1s 2.0 2.5 x 0 
Figure 7. The trapezoidal method for different step-sizes. 


While numerical results for h = 0.5 and 1.0 are not very good they do converge to 
the limit y = 1000 as x increases and reasonably approximate the shape of the 
logistic curve. 


So what about the predictor-corrector method? It is a combination of the Euler 
and trapezoidal methods. Unfortunately, as indicated in the graphs in Figure 8, 
this method also behaves badly for larger values of h. 
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04 08 12 1.6 20% 


Figure 8. The predictor-corrector method for different step-sizes. 


The numerical results follow the shape of the logistic curve for h = 0.1 and 

h = 0.15, but for h = 0.2 the numerical approximations are converging extremely 
slowly towards the limit value of 1000 while for h = 0.25 the numerical results 
converge rapidly to the wrong limit. 


Exercise 3 


Use the predictor-corrector method (Euler-trapezoidal) with h = 0.1 to determine 
approximate solutions for y at x = 0.1 and x = 0.2 for the differential equation 


y' =logy + x with y(0) = 1. 
Quote your results to 3 decimal places. 


Exercise 4 


Use the predictor-corrector method (Euler-trapezoidal) with h = 0.2 to determine 
approximate solutions for y at x = 0.2 and x = 0.4 for the differential equation 


y=siny with y(0) = 1. 

Quote your results to 3 decimal places. 

Exercise 5 

Eliminate Y*,, and Y¥*{, from the recurrence relations for the predictor-corrector method 
Yh. =Y, AY, where Y;, = m(x,, Y,) 


h 
and Yor = Vet 5(¥r + Yess) where ¥F41 = m(Xp+1,¥F41) 


and hence verify that this method is an explicit method. 
[Solutions on pp. 52-53] 
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Summary of Section 2 


1. The television programme looks at several methods of obtaining numerical 
solutions to the logistic equation 


y= 10y(1 = of with y(0) = 100 


including Euler’s method and the trapezoidal method, 

2. A predictor-corrector method uses an explicit method to obtain a crude 
approximation to the solution of a differential equation at x,,, and then uses a 
more accurate implicit method to obtain an improved approximation at x, +1. 

3. The predictor-corrector method based on Euler’s method and the trapezoidal 
method overcomes the difficulty in using the trapezoidal method to obtain 
numerical solutions for non-linear differential equations. The procedure for 
applying this method is given on page 19. The principal term in the local 
truncation error for this method involves h*. 

4. For Euler's method and the predictor-corrector method we need to choose the 
step-size, h, carefully in order to obtain reasonable results. 


3 The analysis of numerical methods 


From the empirical results of Sections | and 2 we have seen that 


(i) _ it is possible to approximate differential equations using a variety of different 
methods, each of which leads to a recurrence relation; 


(ii) the smaller the step-size, h, the better the approximation; 
(iii) larger values of h may give nonsense. 


In this section we want to look at properties of methods theoretically so that you 
can evaluate any new method you might meet. The analysis will also help you to 
choose the step-size required to give reasonable results. 


We shall introduce three new concepts: 


(1) consistency, which is related to (i) above, 
(2) convergence, which is related to (ii) above, 
(3) stability, which is related to (iii) above. 


For simplicity we are only going to consider one-step methods in which the 
differential equation is approximated by a first-order recurrence relation. All the 
methods you have met so far have been one-step methods, but in Section 4 you 
will see a two-step method in which the recurrence relation is of second order (i.e. 
the formula for Y,,., involves Y,—, as well as Y,), The analysis can be extended to 
multi-step methods but we do not have time to do so in this unit. 


3.1 Consistency 


One of the fundamental requirements of any numerical method that we use to find 
approximate solutions to a differential equation is that the recurrence relation to 
be used should be a reasonable approximation to the differential equation itself. In 
other words the recurrence relation must be consistent with the differential 
equation. If this requirement is not satisfied we cannot hope that the solution 
obtained using the recurrence relation is a reasonable approximation to the 
solution of the differential equation. 


As in the definition of the local truncation error, consistency is concerned with the 
local behaviour of the method. We assume that Y, lies on some solution curve so 
that Y, = y,. For implicit methods we also assume that any other values on the 
right-hand side of the recurrence relation are correct. 
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Any one-step method, including the predictor-corrector method of Section 2, must 
involve a recurrence relation of the form 


Year = Ye + holx,, Yes Yer1sh) (dy 
where the function @ depends on the method used. For example the Taylor series 
method of order 2 uses the recurrence relation 


he 
Yer = th¥r + ¥r 


hyn 
=yY, aly + 3¥3) 
and so, for this method, we have 
h 
P(r» Yrs Vrrieh)= Ye +5 Vr. 
Using the condition that ¥, = y, and that the other right-hand side values are 
correct in Equation (1), we have 
Yorn =Vr + hbX rs Vr Vre ish) 
Rearranging this equation gives 


Yi+41—Jr 
ae = (Xr Ve Ver ish). (2) 


The left-hand side of Equation (2) represents the slope of the line joining the 
points (x,,y,) and (x,+1, ¥,41). As A gets smaller and smaller we require that 
Equation (2) should be a better and better approximation to the differential 
equation 

y= m(x, y) (3) 
at x,. This condition holds when h = 0 provided that 

PXr Vr Vrs 0) = m(X,5 Vr). 


(Since y,.1 = y(x, +h) we have y,,, = y, when h = 0.) We thus have a fairly 
simple test for consistency. 


A one-step method is said to be consistent with the differential The definition of consistency 

equation (3) if the recurrence relation can be expressed in the can alternatively be stated as 

form of Equation (1) with lim Y4 i Lie tee 
ho 


P(X VrsVrr0) = m(x,,¥,). 


Example 1 
For the Taylor series method of order 2 we have 


h 
P(X Ves Yetaeh) = Ve + pid 


so that putting h = 0 in this equation gives 
P(Xr Vr Vrr0) = Yr 
= m(X,,¥,) 
as required. 


Hence the method is consistent with the differential equation (3). 
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Exercise 1 
Which of the following methods are consistent with the differential equation 


y =m(x,y) 2 


h 
(i) Yer = Ye + (Vr + Via) 
es h 
Gi) Yrrr = ¥% + GGYs — Yous). 


Exercise 2 
What condition must be imposed on the values of a and b for the method given by 


Your = ¥, + ha, + bY, 41) 
to be consistent? 
[Solutions on p. 53] 


3.2 Convergence of one-step methods 


In Exercises 6 and 9 of Section 1 you were asked to compare the results of using a 
method with two different step-sizes. In Exercise 6 you saw that halving the step- 
size for the Taylor series method of order 2 reduced the errors by a factor of 
approximately 4. In Exercise 9 you obtained a similar result for the trapezoidal 
method, 


What would have happened if we had halved the step-size again? Well, 
commonsense tells us we could expect a further reduction in the error by another 
factor of 4 each time. Thus if we continue halving the step-size (and doubling the 
number of calculations) we would theoretically obtain answers to any accuracy we 
require. (In practice there is a limit to the accuracy that we can achieve because of 
the build-up of rounding errors in the calculations.) Suppose we assume that all 
calculations are carried out exactly, Then our empirical evidence from Exercises 6 
and 9 in Section 1 indicates that the solutions converge in each case. That is, if we 
take some particular point x* (say x* = 0.2) and compute the approximation to y 
at that point using smaller and smaller step-sizes, then these approximations get 
closer and closer to the true solution at that point. Figure 1 illustrates how the 
numerical solution might converge to the true solution. 


ere 
qo h = 0.025 


y= 0.05 
. 
“——h=0.1 
T T T = 
0.1 x =0.2 03 x 


Figure 1. Graph showing how the numerical solution could converge to 
the solution of the differential equation. 


For a given numerical method let us look at the error at a fixed value of x as h 
changes. 


(i) Suppose we take a fixed value of x, say x*. 
(ii) As h decreases so the number of steps, N*, from xp to x* increases such that 
x* = Xo + N*h. 
The global error at x* is the error 
Yue — y(x*). 


We are now in a position to define convergence. 
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1. You are given a problem of the form 
yp’ =m(x,y) with (xo) given. 


2. You wish to use a one-step numerical method to obtain 
approximations for xy < x <b. 


3, The method is convergent for this problem if, for all points x* such 
that x9 < x* < b, we have 


lim Yy+ = y(x*), where N*h = x* — xo 
hoo 


ie, the global error at x* tends to zero as h tends to zero. 


Convergence is difficult to prove directly, except for very simple differential 
equations. However we can state a rather simple but very useful theorem for one- 
step methods which makes any direct proof of convergence unnecessary. 


Theorem 1 
For a given problem 
y' =m(x,y) 


and a value for y(xo), where m is well-behaved, a one-step method 
is convergent if and only if it is consistent. Furthermore if the 
principal term in the local truncation error involves h?*', for 
some integer p, then the global error at x*, for small h, is of the 
form 


Yy- — y(x*) ~ Ch? 
where C does not depend on h. 


This is an important theorem because it gives us a way of evaluating our methods. 
For example, the trapezoidal method has a local truncation error whose principal 
term involves h*, Our theorem tells us that the global error is of the form 


Yue — y(x*) = Ch?, 


While this in itself does not help to put a bound on the error (since we do not 
know C) it indicates that if we were to halve the step-size h then we could expect 
the error in Yy- to decrease by a factor of four, even though we've had to do twice 
as many calculations to compute Yys. This is precisely the result we deduced in 
Exercises 6 and 9 in Section 1. 


The above results indicate that we could classify methods in terms of the principal 
term in their local truncation errors. A preliminary rank ordering of the methods 
we have met so far, indicating their practical limitations, would be as follows. 


Well-behaved here means that 
we can differentiate m as often 
as we wish. 
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Principal term 


involves 
i$ Taylor series method 
of order 3 
Drawbacks: 
second and third 
derivatives 
required. 
ns Predictor-corrector Taylor series method Trapezoidal method 
(Euler-trapezoidal) of order 2 
Drawbacks: Drawbacks: Drawbacks: 
none second derivative for linear 
required, equations only 


Euler’s method 


Drawbacks: 
none, 


This diagram makes the point that if high accuracy is required then the best 
method is likely to be a Taylor series method provided the derivatives can be 
obtained. If less accuracy is required or the derivatives are difficult to compute 
then the predictor-corrector method may be preferred. 


Exercise 3 

Show that the trapezoidal method is consistent and hence convergent, How does the global 
error behave at a point x* for small values of h? If the step-size were halved, by what factor 
would the global error at x* be reduced? 


Exercise 4 
Use the results of Exercise 5 in Section 2 to show that the predictor-corrector method 
(Euler-trapezoidal) is consistent. 


[Solutions on p. 53] 


3.3 Stability 


There is still an outstanding practical issue in the analysis of one-step methods in 
that we do not yet know, given a particular differential equation, what value for 
the step-size, h, to use. The convergence theorem has only told us that, as h tends 
to zero, the numerical solution tends to the true solution—and only then if we 
carry out our computations accurately. For a chosen step-size we need to know 
whether we can expect reasonable results. 


In the television programme we saw examples of how numerical solutions can go 
haywire if the chosen step-size is too large. In each case problems arose because 
errors, incurred in approximating the differential equation, swamped the solution. 
In this subsection we are going to investigate a fairly simple test equation to see 
when these difficulties arise and in the next subsection we will apply our results to 
more practical problems. The following exercise illustrates what can happen. 


Exercise 5 
Use Euler's method with h = 0.1 to solve the differential equation 
y= —30y with y(0) = 1. (1) 


Compare your results with the true solution, y = e~ °°. 


[Solution on p. 53] 
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In Unit 1 two types of ill-conditioning were defined: absolute ill-conditioning and 
relative ill-conditioning. In this unit we are only going to discuss absolute ill- 
conditioning, 


First we consider the differential equation itself. Some differential equations are 
absolutely ill-conditioned in the sense that small changes in the data result in 
much larger changes in the solution, Whatever numerical method we use we 
cannot hope to avoid the build-up of errors, although we can hope that those 
errors will not swamp the solution. In order to deal with this case thoroughly we 
would need to consider relative ill-conditioning, which is beyond the scope of this 
unit. On the other hand there are some differential equations which are well- 
conditioned and it is these problems which particularly concern us here. 


When we apply a numerical method to a well-conditioned problem the recurrence 
relation problem can be absolutely well- or ill-conditioned, In order to discuss 
this we need a new definition which applies to numerical methods rather than to 
the problems they are intended to solve. 


A numerical method, applied to a differential equation, is 
absolutely unstable if the resulting recurrence relation problem is 
absolutely ill-conditioned. The method is absolutely stable if the 
recurrence relation problem is absolutely well-conditioned. 


To illustrate these ideas we consider the test problem given by 
y' =ay with y(0) = 1 


where « is a known constant. This equation can be solved analytically (using, for 
example, the separation of variables method) giving the solution 


y=e™. 
Figure 2 shows the solution for various values of «. 


y 
t a=2 


0 05 
Figure 2. Graphs of e** for various values of 
Suppose we perturb the initial condition for this differential equation so that 
j0) =1+e 
where ¢ is small. The new solution is 


= (1+ ele™ = e* + ce™ = y + ce, 


The error in the solution (¢e**) is significantly larger than the error in the data (e) 
whenever e** is significantly greater than one. In this case the differential equation 
is absolutely ill-conditioned, Now for sufficiently large positive x this condition 
holds whenever wis positive. We conclude that the differential equation problem is 
absolutely ill-conditioned if > 0. By a similar argument the problem is absolutely 
well-conditioned if x < 0. 
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So what happens if we try to use a numerical method to solve the test problem, 
Equation (1)? To illustrate the approach we look at Euler’s method which uses the 
recurrence relation 


Yiii=Y¥, + AY,. 
From the differential equation, Y; = «Y,, so we have 
Yp41 = Y, + hay, 
=(1+ha)¥,. 
We have a recurrence relation of the form 
Yi41=aY, (2) 
where a = 1 + ha. 


The recurrence relation problem is thus to generate a sequence of numbers 
Y,, Yo, Ys... using the recurrence relation (2) with Yo = 1. 


What we now ask is ‘when is this problem absolutely ill-conditioned?” 
From Unit 1 we know that absolute errors grow whenever 

la| > 1. 
(For example if we put Yo = 1 + ¢ then 

Y, = (1 + 8)a" =a" + 6a" = ¥, + ca” 


and the absolute error, g, in Yo has induced an error in ¥, of ea". This error is 
larger than the error in the data if |a| > 1.) 


There are two cases to consider: 

Case (i) a >0 

Since the step-size h is positive we automatically have 
a=1+he>1 


and the recurrence relation problem is absolutely ill-conditioned. This is not 
surprising as the original differential equation problem is absolutely ill- 
conditioned. However, this ill-conditioning may not matter as the growth of the 
relative errors may be tolerable. 

Case (ii) « <0 


This is the more interesting case since we know that the differential equation is 
absolutely well-conditioned, 


For negative « we have 
a=1+ha<1. 


However this does not imply that the method is absolutely stable. For example, 
suppose that h and « are such that ha = —3. Then 


a=1+he=1-3=~2 
and, since |a| > 1, the recurrence relation problem is absolutely ill-conditioned. 
In general the condition for absolute ill-conditioning is that the modulus 

|1 + hol > 1 


so that whenever 1 +ha < —1 we have an absolutely ill-conditioned recurrence 
relation problem. The condition for absolute ill-conditioning is that 


ha < —2, 


Here we have an example in which the original problem is well-conditioned but 
the numerical method gives rise to a problem which is absolutely ill-conditioned if 
ha < —2, 


Figure 3. While absolute 
errors are increasing, relative 
errors may not increase 
significantly. 
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We conclude from this test problem that Euler's method is absolutely stable if h 
and « are such that 


—2<ha <0, (3) 


The set of all values of ha which satisfy this inequality is denoted by (—2,0) and is 
called the interval of absolute stability for Euler’s method. 


We can now use this analysis to understand what went wrong in Exercise 5. There 
we tried to solve the differential equation 


y' = —30y with y(0) =1 
using Euler's method with a step-size h = 0.1, The theory tells us that the method 
is absolutely stable for this problem only if 

—2< -—30h <0, 


ie. h < 1/15. While the theory does not tell us that we will get accurate answers, 
for h = 0.05 say, it does warn us that it is impossible to even get sensible answers 
for the step-size we tried, h = 0.1. 


The tape session in the next section will help you to prove some of the stability 
results for other methods. 


Exercise 6 
For what step-sizes is Euler’s method absolutely stable for the differential equation 


y’ = —100y with yO) =17 
(Solution on p. 53] 


34 Stability in practice 


I mentioned in the last subsection that the absolute stability results we obtained 
for the equation 


y =oy with y(0) = 1 


could be extended to more difficult problems. Indeed the extension to the general 
linear differential equation 


vy! = I(x)y + k(x) 
is quite straightforward. Euler's method gives the recurrence relation 
Yor = ¥,+hY¥, 
= ¥, + h[l(x,)¥, + k(x,)] 
= [1 + Al(x,)]Y, + hk(x,). 
This non-constant-coefficient recurrence relation is absolutely stable (see Unit 1, 
Subsection 3.2) if 
|1 + Al(x,)| < 1, 
ie. —2 <hi(x,) <0 for all the values of r. 
We would thus have to check, before implementing the method, that this 
relationship holds for all the values of x, we are going to use. However, if I(x,) > 0 


we can still use Euler's method with some confidence provided that we know that 
the solution is increasing in magnitude. 


The treatment of non-linear equations is not quite so straightforward. We first 
consider equations of the form 


y= m(y). 
For example 
y' =10y/1- 2, with y(0) = 100, 
1000, 


(This is the problem considered in the television programme.) 


In M101, Block 1, the 
notation Ja,b[ was used for 
the open interval a < x <b. 
Here (—2,0) means the same 
as ]—2,0[, 


This is the same as Equation 
(3) with « replaced by I(x,). 
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If we apply Euler’s method to this type of problem we obtain the recurrence 
relation 


Yiu. = ¥, + hm(Y,). (4) 


To see whether errors grow, we look at what happens to Y,,., if there is an error, 
e, in Y,, so that the computed value for Y, is 


Y=¥,+a 
Then the computed value for Y,,; is 
Fri = + hm(%) 
Y, +8 + hm(Y, + €). (5) 


We can approximate m(y) near Y, by its Taylor polynomial of degree one as 
m(¥, +) = m(¥,) + oy) 
dy 
Hence, using this approximation in Equation (5), we have 
Yur =¥,+et+ hm(y,) + ne™(y,), 
dy 
Using Equation (4) we can reduce this to 
dm 
Bar = You (1 + nin). 
y 
So the error in the computation of ¥,,, is approximately 
dm 
(i + wit). 
This means that the absolute error, «, in Y, has been magnified by a factor 


i nmr) in the computation of ¥.,,. The recurrence relation problem for 


Euler’s method is thus absolutely ill-conditioned if 


14nd (y,) eat 


dy 
Hence Euler's method is absolutely stable if 


dm 
Soh Oreo! (6) 


If this condition is true for all values of Y, in the calculation then Euler's method 
will be absolutely stable for the differential equation y’ = m(y). 


Example 2 
For the logistic equation considered in the television programme we have 


ia Biv! 
m(y) = 101 7000)" 


Differentiating with respect to y gives 


Now y takes values between 100 and 1000 so that = takes values between 8 and 
—10. 

Our conclusions would be as follows: 

(i) For 100 < Y, < 500, amy, is positive and so Euler's method is absolutely 


unstable for any value of h. However, since the solution is increasing rapidly in 
magnitude, this instability may not cause significant errors in the solution. 


This is the same as Equation 
(3) with a replaced by to 
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(ii) For 500 < ¥, < 1000, 3 (¥,) is negative and can be as small as — 10. 


Equation (6) is satisfied with = = —10 only if h < 0.2, Thus we would only 
contemplate using Euler’s method with step-size h < 0.2. 


The above analysis agrees precisely with the results we obtained in the television 
programme. The instability is most serious for values of Y, near 1000. In Figure 6 
of Section 2 we saw that, for values of h less than 0.2, the sequences generated by 
the recurrence relation converged to 1000 while, for h > 0.2, the sequences 
generated did not converge to 1000. This was caused by the absolute instability of 
Euler’s method for h > 0.2. 


What we have found so far, for non-linear equations of the form 

y' =m(y), (7) 
is that the stability condition on the values for h in Euler’s method can be 
obtained by substituting the most negative value ae for « in the inequality 
—2<ha <0. 


It can be shown that, for any numerical method applied to the differential 
equation (7), the stability condition on the values of h can be obtained by 


substituting the most negative value of o for in the inequality which determines 
the interval of absolute stability for that method. 
Finally we consider the general equation For example 
y= m(x,y) aise 
with yo given. with y(0) = 1 for0< x < 10. 
If we apply Euler’s method to this problem we obtain the recurrence relation 
Yin. = Y, + hm(x,, ¥;). 
This can be written as 
Yi+1 = ¥, + hm,(Y,) (8) 
where m, is the function defined by hte eles aa 3 
m,(y) = m(x,,y). function of y only. 


The treatment of Equation (8) is now identical to that given to Equation (4) and 
we can deduce that Euler's method is stable if 


dm, 
rs AY '0' 


for all values of r, 


-2<h 


Hence, when we consider the stability of a numerical method applied to the 
differential equation 


y' = m(x,y), 
we replace « in the inequality which determines the interval of absolute stability by 
dm, 
dy 
previous examples are special cases of this general first-order differential equation.) 


the most negative value of (Y,) where m,(y) = m(x,,y). (Note that all the 


Example 3 
For the equation 


y! = —xy? with y(0) = 1 
for 0 < x < 10 we have 


m,(y) = —x,y? 
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dm, 
dy 


From, for example, a sketch of the direction field for this equation we could show 
that y lies between 0 and 1. Thus, since 0 < x < 10, the most negative value of 
—2x,Y, is —20 and Euler’s method is certainly absolutely stable if 


—2< -20h<0 
ie. h <O.1. 


and 


(¥,) = —2x-¥,. 


From a practical viewpoint, faced with a differential equation, a simple procedure 
for obtaining numerical solutions to a required accuracy is to 


(i) determine bounds, if any, on the values of h to ensure the absolute stability of 
the method; 


(ii) solve the differential equation problem using the numerical method with two 
different step-sizes, h. A comparison of the results at specific x values will 
normally indicate the accuracy. If the discrepancies are larger than desired, a 
smaller step-size should be used and the numerical results compared. 


Exercise 7 
Euler’s method is to be used to solve the differential equation 


y= —xy+ I with y(0) = 1 


+x? 
for0 <x < 10. 
Determine a bound on the values of h which will guarantee that Euler's method is 
absolutely stable. 


Exercise 8 
Euler’s method is to be used to solve the logistic equation 


with the initial condition y(0) = 2000. 


The theory tells us that the solution should tend to the stable population of 1000. 
Determine a bound on the values of h which will guarantee that Euler's method is 
absolutely stable, 


Exercise 9 
Euler’s method is to be used to solve the differential equation 


with y(0) = 1. 


It is known that y takes values between | and 5 for 0 < x < 2. Determine a step-size h for 
which Euler's method will be absolutely stable, 


[Solutions on p. 54) 


Summary of Section 3 
1. A one-step numerical method, which uses the recurrence relation 
Year = Ye + h(x, Yes Yea ish) 
is consistent with the differential equation 


yi =m(x,y) 


PUXr Yes Vr» O) = mUXp, Yr) 
2. The global error at x* is the error 
Yw«— ylx*) 


where y(x*) is the true solution at x* and Yy- is computed using a recurrence 
relation N* times, where 


x* = Xo + N*h, 


If we knew the soluti 
could get a less restrictive 
bound on h but we have 
managed to obtain a useful 
bound on h even without 
knowing the solution. 
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3. A one-step method is said to be convergent for the problem 
y’ =m(x,y) with y(x9) given, 


if the global error at x* tends to zero as the step-size, h, tends to zero, for all 
points x* in the interval, [xo, 6], in which approximations are required; 


ie. lim ¥ye = y(x*), 
h=0 
where N*h = x* — xo. 


4. A one-step method applied to a well-behaved differential equation is convergent 
if it is consistent (see Theorem 1). If the local truncation error involves h’*' 
then the global error at x*, for small h, is of the form 

Yy+— y(x*) = Ch? 
where C does not depend on h. 

5. A one-step method, applied to a differential equation, is absolutely unstable if 

the resulting recurrence relation is absolutely ill-conditioned. The method is 


absolutely stable if the recurrence relation problem is absolutely well- 
conditioned. 


6. The interval of absolute stability (a,b) is the set of all values of ha for which the 
method is absolutely stable. 


7, For Euler's method we have the following intervals of absolute stability 


differential equation | interval of absolute stability 


y= ay —2<ha<0 
y' = I(x)y + k(x) —2<hl(x,)<0 for all r 
y= m(y) -2 < iH) <0 for all r 
dm, 
y! = m(x, y) RST ore iO Bac where m,(y) = m(x,, y) 


For other methods it is sufficient to determine the interval of absolute stability for 
the test problem 


y =ay with y(0) =1 


and then to deduce the appropriate conditions for other differential equations as 
above. 


4 Exercises on stability and Simpson’s 
method 


4.1 Exercises on stability (Tape Subsection) 


In Subsections 3.3 and 3.4 we looked at the absolute stability conditions for 
Euler’s method. In the first part of the tape session we see how the intervals of 
absolute stability are determined for the trapezoidal method and the predictor- 
corrector method. 


In the second part of the tape we apply these results to the logistic equation of 
Section 2 to confirm the empirical evidence found in the television programme. 


Start the tape when you are ready 
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y=xy y(0) =1 


Recurrence relation for the trapezoidal method 


h ’ ' 
Yee 2 Ce +Y/44) 


Graph of 
siete 


la] <1 when ha < (2a) 


The interval of absolute stability for the trapezoidal method is (—00, (=) 
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lal<lwhen]  |<he<[ | (4a) 


The interval of absolute stability for the predictor - corrector 


method (Euler -trapezoidal) is Ppsare a 
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The logistic equation 


If 500 < y < 1000 the minimum value of a is eid 


The trapezoidal method 


The interval of absolute stability is (—<0, 0) 
.. the method is absolutely stable if he < 0. 


For the logistic equation, with 500 <y < 1000, what 


values of h make the trapezoidal method absolutely stable ? 


re 


Lz) The predictor ~corrector method (Euler ~ trapezoidal) 
The interval of absolute stability is (-2,0) 
The method is absolutely stable if -2 < he <0 
For the logistic equation the method is absolutely stable 
for 500 <y < 1000 if 


-2< Was <0 
-2<-10h 
ieh <0-2 
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Exercise 1 
Find the interval of absolute stability for the Taylor series method of order 2, applied to the 
test equation y’ = ay with yO) = 1, given by 

he 


Yer = Hrhh +> kr 


where Y7 =aY, =27Y,, 
For what values of h is this method absolutely stable when applied to the logistic equation 


y= 101 = ita with y(0) = 500? 


[Solution on p. 54) 


4.2 Simpson’s method 


For simplicity we have only analysed one-step methods so far. Most of the 
methods used in practice are many-step methods. To give you some idea of how 
they work we are going to consider a two-step method: Simpson's method. Like 
the trapezoidal method Simpson’s method is an integration method. Because of its 
high accuracy Simpson’s method was very widely used, before the advent of 
computers, to solve differential equations. It is possible to show that the local 
truncation error for Simpson’s method involves h* so that the global error using 
this method is likely to be much smaller than for most of the methods we've 
already discussed. 


Simpson's method approximates an integral over two steps as 
re h 
{| S(x)dx = 5 (f(%r—1) + MC) +S %r41)). (1) 
rot 
We are going to use this formula to solve the differential equation 
y= m(x,y)- (2) 
If we integrate both sides of Equation (2) over the interval [x,_,,x,41] we have 
eet ee 
f yidy= f m(x,y) dx. 
emt 


enn 
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The left-hand side is just y,,, — y,—1 So that, using Simpson’s method (1) for the 
integral on the right-hand side, we have 


Yeti Vena -{ m(x,y) dx 


h 
5 (Im(X,—15e—1) + Amare) + M(Xp+1sVr+1))- 
This leads us to define the recurrence relation for Simpson's method as 
h 
Yor = Yr-1 +3(Ki-1 +4Y, + Yi4s) 


where Y’, = m(x,, Y,). 


This is a second-order recurrence relation (because it involves terms at x, and 
x,~1) and it is also an implicit method (because it involves terms in Y,,, on the 
right-hand side). This causes two difficulties: 


(i) Second-order recurrence relations require two initial values, Yo and Y,, in 
order to generate the sequence Y>, Ys, Y,... . However, we only have one 
given initial condition, Yo, so that somehow we have to generate a value for Y, 
in order to use Simpson's method. This is done by using another method, such 
as a Taylor series method, to generate Y, before switching to Simpson's 
method. 


(ii) Simpson’s method cannot be easily used on its own to solve non-linear 
equations because it is implicit (cf. the trapezoidal method). For this reason 
Simpson’s method is more commonly used as a corrector in a predictor- 
corrector method. 


By a process of algebraic manipulation it is possible to obtain a linear recurrence 
relation for Simpson’s method applied to the linear differential equation 


y= I(x)y + k(x) 
with y(xo) given. 


I will not present those manipulations here but the process is similar to the 
procedure for the trapezoidal method in Section 1. The result is given in the 
following procedure box. 


r Simpson’s method for linear differential equations 
1. You are given a linear differential equation 
y= Mex)y + kx) 
and a value for y(xo). 


2. To apply Simpson's method, choose a step-size h and calculate 
Y), Ys,... from the second-order recurrence relation 


_ ({_Anl, 3+hl,-1 h 
Yai = (a), + Gta). +(e. + 4k, + kes) (3) 


where Yo = »(Xo), ¥; has been previously obtained using some 
other numerical method and x, = Xo + rh ly 


3. Y, is an approximation to y,. 


In the following example we use this procedure to illustrate the technique, 


Example 1 
Use Simpson’s method to solve the differential equation 
y=axty with y(0) = 0 (4) 
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using a step-size h = 0.1 for0 <x <1. Initially generate the value of Y, at x = 0.1 
using a Taylor series method of appropriate order. 


Solution 

Since Simpson’s method has a local truncation error which involves h° it is 
worthwhile computing Y, very accurately. The Taylor series method of order n 
uses the recurrence relation 

W hn 

J yr 4 eg Dy 
5 Yopoy 7 yo 


where the derivatives may be obtained using the differential equation (4). We have 


Year = Y¥, +hY; + 


Yi= x, + ¥,; Yr=1+Yj, yr=Y; and so on. 
With h = 0.1 we have 
1? |, , 0.1)? 
Dieta 


where % =0, Yj =xo + % =0, Yo=1+Y¥o=1, Yo'=Y$=1 andsoon, 


0.1)" 
Y=%+01% + Ye" te OT yp 


os @Ay  ©.1)P | Ole uy" 
Hence ¥, =0+0+ 7+ 6 + 57 foee$ 


= 0.005 + 0,00016667 + 0,00000417 + 0.00000008 + ++ 
= 0.00517092. 


The appropriate Taylor series method for eight decimal place accuracy for Y, is of 
order 5, 


With this value for Y, we use the procedure for Simpson's method, We have 
U(x) =1 
and k(x) = x. 
Hence the recurrence relation (3) becomes 
Teste (4) ¥+ (544) Bats (rgles + 48, + Xp) 
With A = 0.1, and noting that 
Xp-1 + Xp4 = 2X, = Irh, 


we have 
4 31 1 
39 y+ 29 Y-1 + 2» 
This second-order recurrence relation can be used to generate the solution using 
the initial conditions Y) = 0 and Y, = 0,00517092. In the following table the true 
solution 

yoe—-x-1 
is also given for comparison. 


Yer = (0.6r). 


r oF ¥, (Simpson’s method) true solution error 

0 0 0 0 0 

1 0.1 0.00517092 0.0051709169 3 x10 
2 0.2 0.02140289 0.02140276 1.3 x 10-7 
3 03 0,04985897 0.04985881 1.6 x 10-7 
4 04 0,09182501 0.09182470 Sie 1057 
5 0.5 0.14872166 014872127 3.9 x 1077 
6 0.6 0.22211938 0.22211880 5.8 x 1077 
7 0.7 0.31375341 0.31375271 7.0 x 10-7 
8 0.8 0.42554187 0.42554093 94x 10-7 
9 09 0.55960425 0.55960311 1.14 x 107° 
10 1.0 0.71828328 0.71828183 1.42 x 10°° 

al 


You may get slightly different 
results on your calculator, 
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All the results generated using Simpson’s method are correct to five decimal 
places. Comparing these results with the results obtained in Example 1 of Section 
1 for the same differential equation it can be seen that the results are much better 
for Simpson’s method than for Euler's method, the trapezoidal method or the 
Taylor series method of order 2. Although quite a lot of work was necessary to 
compute a value for Y;, the results clearly make this effort worthwhile if great 
accuracy is required. 


Exercise 2 
Use Simpson's method to solve the differential equation 
y’ =3y + sinx with y(0) = 0. 


With h = 0.2, compute ¥; and Y, given a value Y, = 0.0255, (¥, was obtained using the 
trapezoidal method with a step-size h = 0.1.) 


Exercise 3 

Show that the use of Simpson's method, with h = 0.1, applied to the differential equation 
y’ = 10y with yO) = 1 

leads to the linear second-order recurrence relation 
b PP) Phe) eres 

By solving the auxiliary equation 
x? =2x4+2 

write down the general solution of the recurrence relation. 

[Solutions on pp. 54-55) 


4.3 The numerical analysis of Simpson’s method 


As you saw in the last subsection, Simpson’s method does give very accurate 
results provided we obtain a good approximation to Y, using some other method. 
If hundreds or thousands of calculations are to be carried out then Simpson's 
method is likely to give much better results with fewer calculations than any of the 
high-order Taylor series methods. 


However, there is one further drawback associated with Simpson’s method which 
we hinted at in Exercise 3: Simpson's method cannot be used to solve certain 
problems. The use of Simpson’s method for a linear equation leads to a second- 
order recurrence relation. From Unit 1 we know that the general solution of a 
second-order recurrence relation contains two arbitrary constants. For example the 
general solution of the recurrence relation in Exercise 3 can be written as 


Y= Ax Aan+ Bx pl" 
where A = 2.7320508 and p = —0.73205081. 


In Unit 2 we saw that the general solution of a first-order recurrence relation 
involves only one arbitrary constant. For example the general solution for the 
differential equation 


y' = 10y 
is y= Ae™, 


This general solution contains only one term while the general solution for 
Simpson's method contains two terms. What is the correspondence between these 
two solutions? The general solution to the differential equation can be written as 


y, = AelO% 
and 

Ven = Actors 

= Ae! +m 


= Ael™ x gioh 


= ely, 
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Hence there is a recurrence relation for generating y,,, from y,. With h = 0.1 we 
have 

Veri = Vr 

= 2,7182818y,. 

The general solution of this recurrence relation is 

Vn = A(2.7182818)", (5) 
Compare this with the solution of the recurrence relation 

Y, = A(2.7320508)" + B(—0.73205081)". (6) 


There is a similarity between the first term in Y, and y,. The second term in ¥,, 
B(—0,70205081)", is called a spurious solution because it is an extra solution 
introduced by the method. It does not correspond to anything in the original 
differential equation. How will this affect the results? 


We will apply Simpson’s method to the differential equation 
y’ = 10y with y(0) = 1. 
Suppose that we start the Simpson calculations by computing ¥, = 2.718 using 


some other numerical method, Using this value and Yj = 1 we can calculate the 
appropriate values for A and B in (6). The particular solution is 


Y, = 0.9959 x (2,7320508)" + 0.0041 x (—0,73205081)". 
The particular solution for the differential equation is obtained by substituting the 
initial condition, yp = 1, into Equation (5) giving 

Vy = (2.7182818)". 


Note that the first term in ¥, is almost identical to the true solution, y,, while the 
spurious solution has a small coefficient (0.0041). As n increases the spurious 
solution gets smaller and smaller since (—0.73205081)" tends to zero and hence its 
effect on the solution diminishes. We can conclude that in this case the spurious 
solution introduced by Simpson's method does not seriously affect the results. 


In Exercise 4 you will see that, when Simpson's method with h = 0.1 is applied to 
the differential equation 


y' = —10y with y(0) = 1, 
the spurious solution does seriously affect the results. Without going into the 
details of the analysis we state that Simpson's method should not be used to solve 
differential equations of the form 

y! = Ixy + k(x) 


if I(x) is negative for any value of x in the interval under consideration. This is a 
generalization of the results obtained from Exercises 3 and 4. 


Exercise 4 

Use Simpson's method to write down the recurrence relation for the differential equation 
y= —10y 

with yO) = 1. 


Using a step-size h = 0.1 show that the general solution of the recurrence relation is 
¥, = A(0.3660254)" + B(—1.3660254)", 

What is the particular solution with Yo = 1 and ¥, = 0.3679 7 

Compare this with the true solution given by 
Yn = (0.36787944)". 

What happens to ¥, and y, as n increases? 


(You will have an opportunity to see the effect of this spurious solution when you use the 
computer package.) 


[Solution on p. 55] 
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Summary of Section 4 


1, The intervals of absolute stability for first-order methods are given in the 
following table. 


interval of 
method absolute stability 
Euler's method (—2, 0) 
trapezoidal method (—o0, 0) 
Taylor series method of order 2 (—2, 0) 
predictor-corrector method (—2, 0) 


2. Simpson’s method can be used to generate approximate solutions to linear first- 
order differential equations using the procedure on page 40. 


3. Simpson's method can give rise to very accurate solutions since the principal 
term in its local truncation error involves h°. However, Simpson's method 
should not be used to solve linear differential equations of the form 


y’ = Ix)y + k(x) 


if I(x) is negative as, in this case, the spurious solution (the extra solution 
introduced by the method) will eventually swamp the solution. 
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5 The computer package NUMSOL 


5.1 Introduction 


This section is a description of the computer package NUMSOL. The package is 
designed to solve first-order differential equations using the methods described in 
this unit. It is intended that you should use this package at Summer School, 
although it will be available if you want to visit your local computer terminal. A 
general description of how to input data to packages is given in Unit 1 Section 6. 
This package is very similar to the package RECREL described in Unit 2. 


To help you to decide which method to use for a particular problem I have 


summarized the properties of each method in the following table: 


principal term 
in the 


method Ite involves 


advantages 


disadvantages 


Euler’s method H 


simple to use; 
explicit method 


not very accurate; 
interval of absolute 
stability is (—2,0) 


Taylor series 


explicit method 


requires computation 


method of of y”. Interval of 
order 2 absolute stability 
is (—2, 0) 
predictor- explicit method interval of absolute 
corrector stability is (—2,0) 


method (Euler- 
trapezoidal) 


trapezoidal ns interval of implicit method 
method absolute stability 
is (— 21,0) 
Taylor series ht explicit method: interval of absolute 
method of very accurate stability is 
order 3 (—2.51, 0)* Requires 
the computation of 
y” and y". 
Simpson’s method ns very accurate implicit method; 


for some problems 


should not be used 
for y’ = I(x)y + k(x) 


if I(x) is negative 


* not given in the text earlier. 


To illustrate how the package works we begin with an example showing how the 
logistic equation can be solved at the terminal. 


Example 1 
Solve the logistic equation 


a pee 
= = 10y(1 an 


with y(0) = 100 


on the interval [0,2] using the predictor-corrector method with h = 0.1. 
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Solution . 
The following terminal dialogue shows how I input this problem, after logging in 
and obtaining the course library program NUMSOL. Each of my responses is 
underlined to indicate that this is what I typed, The other characters in the 
dialogue were output by the computer. 


Option 10 specifies the differential 


OPTION? 10 equation to be solved 


Y' = ?10*Y#(1 —Y/1000, 
OPTION? 23 


PREDICTOR-CORRECTOR METHOD (EULER-TRAPEZOIDAL) 
OPTION? 30 


INITIAL CONDITION Y(X(0)) =? 100 
OPTION? 32 


These two options specify the 
range of x values 


INITIAL X VALUE X(0) = ?0 
OPTION? 33 
FINAL X VALUE X(N) =?2 


OPTION? 34 


STEP-SIZE H = 70.1 


OPTION? SOLVE Having input all the data, you can solve the problem 


5.2. The options for NUMSOL 

Before using the package at a terminal you should plan which options.you are 
going to use to solve your problems. A description of each option available for the 
package is given below. 

Command options 


OPTIONS — this tells the program to print a list of the available options. 


SOLVE —this tells the program to run the problem after all the data has been 
input, If some data is missing an error message will be printed, 

HELP —this tells the program that you are stuck and need advice. 

LIST —this tells the program to print the problem and data previously input 


so that you can check the problem you are solving. 
STOP —this is the only way to stop the program, 


ANSWER —this gives the solution to the computing exercises so that you can 
check that your answer is correct. To the question 


EXERCISE? 


you respond with the number of the exercise you are interested in. 


Problem options 
10 —Input m(x, y) 
To the question 
v=? 


: ; Valid input expressions are 
you respond with a valid input expression for m(x, y). described in Unit 1. 
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For linear equations you must input 
U(x) x y + k(x) 
in that order so that the package can recognize the linearity. 


—Order and derivatives for Taylor series method 
To the question 
ORDER =? 


you respond with the order (<5) of the Taylor series method. When 
you have specified the order the program will ask you for each of 
the higher derivatives that the method requires. 


For example if you want to use a Taylor series method of order 2 
for the differential equation 


y=xyt3 


you would respond to ORDER =? with 2. The package would then 
ask you for the second derivative 


Y=? 
to which you would respond Y + X*Y' 


(The differential equation would have been input using Option 10 so 
Y’ is already known.) 


12 —Change order of Taylor series method 
If you wish to change the order of the Taylor series method without 
having to input previous derivatives again, you would use this 
option. To the question 

ORDER =? 
you would respond as in Option 11 and if further derivatives are 
required you will be asked for them. 
For example if you use the same differential equation as in the 
example in Option 11 but wish to change to the Taylor series 
method of order 3, to the question ORDER = ? you would respond 
with 3. The package already has Y’ (Option 10) and Y” (Option 11) 
and so only Y” is missing. The package will ask for this and to the 
question 
Y=? 

you would respond 2+ Y' + X#Y"’ 

Method options 

20 —Euler’s method 
This uses a first-order recurrence relation to solve the differential 
equation. The initial y value is input using Option 30. 

21 —Taylor series method 
The order of the Taylor series method and any further derivatives 
required are specified using Option 11 and the order can be changed 
using Option 12. The method uses a first-order recurrence relation to 
solve the differential equation. 

22 —Trapezoidal method 


This method can only be used to solve linear differential equations 
using a first-order recurrence relation. 


Note that only single quotes, 
repeated as necessary, are 
used for derivatives e.g. Y'" not 
Mies 
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23 —Predictor-corrector method 
This method uses Euler’s method to predict and the trapezoidal 
method to correct. Only one initial y value needs to be specified 
(Option 30). 
24 —Simpson’s method 
This method can only be used to solve linear differential equations. 
The method involves the use of a second-order recurrence relation 
for which two initial y values are required (Options 30 and 31). Y; 
may be computed using a one-step method with small step-size if it 
is not known. 
Data options 
30 —Initial condition Yo 
This specifies the initial y value. 
31 —Initial condition Y, 
This specifies the second initial condition for Simpson's method. 
32 —Xo 
—This is used to input the initial x value. 
33 —Xn 
This is used to input the final x value. 
34 —The step-size h 
This parameter specifies the step-size for the numerical method, 
Print options 
40 —Outline printout 


If you only want selected output you should use this option to 
specify k such that every kth term in the sequence is printed. 


41 —Full printout (default) A default option is the one the 
; . ; program will use if no other is 
All terms in the sequence will be printed. specified. 
42 —Print solution only 


Only the last term in the sequence is printed. (Intermediate terms 
will not be printed.) 


5.3 Computer exercises 


The following exercises will require the use of the computer package. Before you 
run the exercises at a terminal you should plan how you are going to tackle the 
problems you wish to solve. Judicious preliminary work will enable you to make 
the best use of your limited time at the terminal. You may not have time to solve 
all these problems. 

Exercise 1: This exercise will reproduce the results of the television programme. 


Use the package to solve the logistic equation 
y & 
y’ = 10y{1 — with y(0) = 100 
4 »(1 1000) BOSS 


for0<x <1, 
by each of the following methods 


(i) Euler's method, 
(ii) predictor-corrector method, 
(iii) Taylor series method of order 2. 
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For each method use a step-size h = 0.1. 


If you have time, look at what happens for larger step-sizes. 


Exercise 2: To investigate stability 
The differential equation 
. 1 4 
Y= wt with y(0) = 1 
is to be solved on the interval [0, 10] using each of the following methods 


(i) Euler's method 
(ii) trapezoidal method 
(iii) Simpson's method 


with step-sizes h = 1, 0.5 and 0.1, 
(To reduce output, print the results only at the integer values of x.) 
Do the results agree with the stability analysis in Sections 3 and 4? 


Exercise 3: To investigate accuracy 
What methods can be used to solve the differential equation 
y=-e with y(0) =1 
‘on the interval [0,1]? 
Compare the numerical results obtained using Taylor series methods of order 2 and order 3 
and a step-size of h = 0.1, To get more accurate results use h = 0.05. 
Hint: fe = yl’ = —e, 


Is it more accurate, for this problem, to use a step-size h = 0.05 with a Taylor series method 
of order 2 or to use a step-size h = 0.1 with a Taylor series method or order 3? 
Exercise 4: To investigate stability 
Experiment with the differential equation 
y’ = —30y with y(0) = 1 


for 0 < x <1, to test the stability of the methods in the package for different values of h, 
The stability of Euler’s method for this equation is discussed in Section 3, 


Exercise 5: To investigate stability 

Use Simpson's method with h = 0.1 to solve the differential equation 
y= —10y with y(0) = 1 

for0<x<2. 

Confirm the results of Exercise 4 in Section 4. 

(Take Y, = 0.3679 and then ¥, = 0.3660.) 


[Solutions to the computer exercises are not given in the unit. A check solution can be 
obtained at the terminal using the option ANSWER.) 


6 End of unit exercises 


Exercise 1 


Compute the first 3 terms in the sequence if the trapezoidal method is to be used with a 
step-size h = 0.1 to solve the differential equation 


y=y+2x with y(0) = 1. 
Exercise 2 
The predictor-corrector method is to be used with h = 0.1 to solve the equation 
y= with y(0) = 1. 


Write down a single recurrence relation which can be used to generate the sequence, What 
is the particular solution for this recurrence relation? 


The valid expression for e” is 
EXP(Y) 


Exercise 3 
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h a | 
Is the recurrence relation Y,,, = Y, +3(4Y, + Y},1) consistent with the differential 


equation 


y =m(x,y)? 


Would you expect this method to give better results than the trapezoidal method for small 


values of h? 


Exercise 4 


For what values of ha is the method in Exercise 3 absolutely stable when applied to the 


problem 
ysay 
with y(0) = 1? 


Exercise 5 

Simpson's method is to be used to solve the linear differential equation 
y' =3y + 2x 

with y(0) = 1. 


Derive the recurrence relation which can be used to generate the solution for a given step- 


size h, 


What other information do you require and how would you get it? 


[Solutions on pp. 55-56] 


Appendix 


Solutions to the exercises in Section 1 
1, The trapezoidal method utilizes the recurrence 
relation 


h 
Yur = Yet 5 le + Vrorh 


For the differential equation 
y =3y +sinx 


we have 
h 
Your = Ye + 5(GY, + sinx,) + BYpe1 + sinxy)]. 


Collecting up the terms gives 


3h 3h\ oh. A 
ti[t-F =» +3) gms, + sinxys1) 


_ (1 + 3h/2)¥, + (h/2\{sinx, + sinx,+1) 


Tent 1—3n72 


[Alternatively we could just use the procedure box to write 
down the recurrence relation directly. ] 


With h = 0.2 we have 


1.3¥, + O.1(sinx, + sinx,.1) 
07 4 


Yor = 


The figures in the following table are quoted to 5 decimal 
places. 


0.2 


04 0.6 


0.02838 | 0.13672 | 0.39020 


0,02460 | 0.12308 | 0.35304 


0.00378 | 0.01364 | 0.03716 


Whilst the results are not particularly accurate, because of 
the rather large step-size, the ‘shape’ of the approximate 
solution is correct, The maximum error is approximately 0.04 
at x = 0.6. 


2. The formula for the trapezoidal method is 
h 
Yi =Ye+ a(t + Vivi 


For the non-linear differential equation 
y' = siny 


we have 
he , 
Yoi=¥+ poke + sin Y,4 1). 
Collecting the terms in Y,,, and Y, we have 


LA ie 
You — 7a ees =¥, +3sin¥,. 
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Putting h = 0.1 gives 
Y,+1 —005sin¥,,, = ¥, + 0.05sinY,. 

With y(0) = Yo = 1, to get Y, we need to solve the equation 
Y, —0.0Ssin Y, = 1.04207355. 


The solution cannot be written down immediately and the 
Newton-Raphson method, as described in Unit 18, Section 2, 
would have to be used to determine Y,. This method 
involves a considerable amount of extra work at each step, 


3. Assuming that Y, lies on a solution curve so that 


Y=) Ye=y, and Ye=yy 


we have, for the Taylor series method of order 2, 
We 
Yreu= Yet hy + yr 


The true solution y,,, can be determined using a Taylor 
series expansion at x, as 
Pi es 
r= US di flr oh Ales oles a) 
Vett Yet lt ay 6? + + 
Hence the local truncation error is 
ns ht 
Voea Yeti oe — agrees 


The principal term in the local truncation error is 


4. Assuming that ¥, lies on a solution curve so that 


Y= Yn Ye=¥n Yr=yr and Ye =y/ 


we have, for the Taylor series method of order 3, 
he elites 
Your = yet hy + or + Eye 


The true solution satisfies 
ne n° nt jlivy hs A 
ors = Yet hye + Ve Eve + aay + Toye 
mel 


Hence the local truncation error is 


nt ns 
Fron —Yrei = — 5) — Tee 


The principal term in the local truncation error is 


WF a) 


i ae 


If the step-size is halved so that h* = h/2 the principal term 
in the local truncation error becomes 


(HY ay = — WAY iy ay 
Tg = Saag ne Se 


‘The local truncation error has been reduced by a factor of 
approximately 16. Note however that about twice as many 
calculations would be needed to obtain an approximation at 
a particular value of x so that we would not expect the 
accumulated error at this value of x to be reduced by a 
factor of 16. 


ht 


5. The Taylor series method of order n is given by 


W i 
Yp41=Y¥,+h¥, +—¥7+---+—YO. 
2 n! 
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Assuming that Y, lies on a solution curve so that 


Y= VY YY = 
we have 
h h" 
va hy tay tore by, 
rt = Yet hy, 2) desert = B 4 


The true solution satisfies 


h? 
Yess = Yet hye + ye 


ra pet 
COREY eee lee 
me’ et 
ie 
Hence the local truncation error is 
yer nea 
eee (ott) nt) 
Bas saa 1) id near” 


6. Here is the table of errors, for h = 0.1, using the Taylor 
series method of order 2 to solve the differential equation 


yvooxt+y 
with y(0) = 1. 
x, y Y— Yr 
0 0 0 0 
0.1 0.005 0.005171 —0.000171 
0.2 0.021025 0.021403 —0,000378 
03 0.049232 0,049859 — 0.000627 


0.05 we have the following table; 


oe Y, Ve 


pp 
0 0 0 

0.1 0.005127 0.005171 

0.2 0.021305 0.021403 

0.3 0.049696 0.049859 


At x = 0.1 the errors are —0.000171 and —0.000044, Hence 
the error has been reduced by a factor of 171/44 = 3.89. 
Similarly at x = 0.2 and 0.3 the errors have been reduced by 
factors of approximately 3.86 and 3.85 respectively. 


We can conclude that, for the Taylor series method of order 
2, halving the step-size reduces the errors at particular values 
of x by a factor of approximately 4. (Note that if we halve 
the step-size the local truncation error is reduced by a factor 
of approximately 8.) 


7. Assuming that the values on the right-hand side are 
correct, we have 


Year =Yr + hyrer- 
The Taylor series expansion for y;,,, is given by 


er 
Voor = Vet hy + Tye +P oo. 


i 
Year = Yet hy + Wye + aye to 


The Taylor series expansion for y,,, is 
3 


nusihgtaees, 
Veri=Vr Vr 2" 6" . 


$2 


Hence 
Ae, 
Yer —Yeos = get ayn to - 


The principal term in the local truncation error is 


rr 
pia 


This term involves h?, as does Euler’s method. 


8. (i) The principal term in the local truncation error for 
the trapezoidal method is 

w 

PD” re 
If we halve the step-size, putting h* = h/2, the principal term 
in the local truncation error becomes 

La PS Li 


TR’ at x T”" . 
Hence the local truncation error will reduce by a factor of 
approximately 8, 


The Taylor series method of order 2 also has a local 
truncation error whose principal term involves h* and so 
halving the step-size would reduce the local truncation error 
by a factor of approximately 8 for this method too. 


(ii) From the results of Exercise 6, for the Taylor series 
method of order 2, we would expect the accumulated errors 
for the trapezoidal method to reduce by a factor of 
approximately 4 if the step-size were halved. 


9. At x =0.1 the error has been reduced by a factor of 
92/23 = 4.00. Similarly at x = 0,2 and x = 0,3 the errors have 
been reduced by factors of 4.00 and 4.02 respectively. 


We conclude that, for this problem, halving the step-size 
reduces the errors at particular values of x by a factor of 
approximately 4 for the trapezoidal method. 


Solutions to the exercises in Section 2 
1, Separating the variables (see Unit 2) we have 
ane re 
y(1 = y/1000) 


To get the left-hand side into the required form we use 
partial fractions to write 


10. (1) 


1 A B 
y(1 = y/1000) y/1000)" 
This leads to the two equations for A and B as 
A=1, 
Foot Bao. 
Hence Equation (1) becomes 


pe ee 
yl =y/1000)  y 1000 — y 
Integrating both sides of Equation (1) gives 


logy — log(1000 — y) = 10x + ¢, 


= 10. 


y 


log| 1000 = y, 


)=10+¢ 
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giving 
a 
1000 — y 
This simplifies to 
- 1000 4e!* 
1 + Aelo* 


= Ael™, 


= 
For the initial condition (0) = 100 we have 


100 = 


1+A° 


1 
Hence A = ry 
The particular solution is 


1000e!° 
Y= 94 er 


1000 
Ge 4 1 
Here is a sketch of this solution. 


iey= 


0 5 1.0 


2. Using the predictor-corrector method to evaluate ¥; and 
Y, we proceed through the four steps of the algorithm: 


j Lo eat 
(i) Evaluate Y= ior, (1 it] 


(using the differential equation), 
(ii) Predict Y¥ = Yq + hYs = 672.52255 


= 2447.79 


(Euler’s method). 


= rors(t = al = 2202.3597 


(iii) Evaluate Y¥ 7000) 


(using the differential equation). 


(iv) Correct Ys = Ya + 4%, + Yt’) = 660.25103 
(trapezoidal method). 


Hence our approximation at x=0.3 is Y; = 660.25103. 
Repeating the steps for Y, we have; 


i ae Fenn 
(i) Evaluate ¥5 = 10%4(1 ita) = 243.1961 
(ii) Predict Yt = Ys + hYs = 884.57064. 


re ve og het2 ls 
(iii) Evaluate Y}' = rors(1 7000) ~ 1021,0542, 


(iv) Correct Ys = Ys + aH + Y3’) = 823.46355. 


Hence our approximation at x = 0.4 is ¥, = 823.46355, 
3. For the differential equation 

y' =logy+x with yO)=1 
and h = 0.1, the predictor-corrector method gives: 
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(i) Bvaluate ¥4 = log Yo + xp = logi +0 =0. 
(ii) Predict Yt = Yo + hYo = 1. 
(iii) Evaluate Y¢' = log Yt + x, = 0.1. 


(iv) Correct ¥; = Yo + ne + Y?’) = 1.005. 


Hence Y, = 1.005 at x, = 0,1. 


(i) Evaluate ¥; = log Y, + x, = 0.10498754. 
(ii) Predict Y$ = ¥, + hY; = 10154988. 
(iii) Evaluate ¥$’ = log Y¥ + x, = 021537988, 


(iv) Correct Y= Y, + 3h + Y$') = 1.0210184, 


Hence Y; = 1,021 at x = 0,2. 
4. For the differential equation 
y’ =siny 


with »(0) = 1 and h = 0.2 the predictor-corrector method 
gives the following approximations at x = 0.2 and x = 0.4; 


(i) Evaluate Y= sin Yo = 0.84147099. 
ii) Predict Yt 1.1682942, 
(iii) Evaluate Yq’ = sin Yt = 0.92008374. 


h 
liv) Correct ¥, = Yo + 5(¥o + Yt") = 1.1761555, 


Hence ¥, = 1.176 at x = 0.2, 


(i) Evaluate Yi, = sin ¥, = 0.92313471, 
(ii) Predict Y¥ = Y, + hY; = 1.3607824, 
(iii) Evaluate ¥3' = sin ¥$ = 0.97802802, 


(iv) Correct Ys = Y, + Ay, + Yf') = 1.3662717. 


Hence Y, = 1.366 at x = 0.4. 


5. Yt. = ¥, +h¥, = Y, + hm(x,, ¥,) using the differential 
equation. 


h 
Yor = Ve + 5(¥e + YR) 


h 
= Ye + 5 Lemp, Ye) + (xe 1 YF 4). 
Substituting for Y#,, in this second equation gives 
Yori =¥,+ t [m(x,, ¥,) + m(x,-, ¥, + hm(x,, ¥,))] 


As this does not have Y;,, on the right-hand side the 
predictor-corrector method is explicit. 


Solutions to the exercises in Section 3 
1. (i) b0%, Yes Yee A) = yet 2Y;41) 
Hence 

(Xr Ver V0) = ior + 2y,) 


= Vr = (Xp Vr) 


Thus the method is consistent with the differential equation. 
1 

(i) P01 Yes Yrs) = GY; — Yrosh 

Hence 


1 
Pr Yer Yer0) = 4 Be — yr) 


ie 
= pir = amr Ir)- 
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Thus the method is not consistent with the differential 
equation. 


2. P(x, Fr, Yury) =aY, + bYy 41. 
Hence 
P(%,,¥e5¥r, 0) = ay; + bys 
=(a + b)m(x,, 


The method is consistent with the differential equation 
y= m(x,y) 
only ifa+b=1, 
If a= b =} we have the trapezoidal method. If a = 1, b= 0 
we have Euler’s method. 


3, From Exercise 2 with a = b = 4 we know that the 
trapezoidal method is consistent. From Theorem 1 the 
method is also convergent. Since the principal term in the 
local truncation error for the trapezoidal method involves h°, 
Theorem | tells us that the global error at x* is of the form 


Yy«— y(x*) = Ch? 
for small values of h, Thus halving the step-size would reduce 
the global error at x* by a factor of approximately 4, 
4, From Exercise 5 in Section 2 we have 


Yui=Y+ A Lm(x,, ¥,) + (xp, ¥, + h(xp, ¥,)) Je 


Hence 
DU% rs Yon eras h) = SLmlx,, ¥,) + m(x,, ¥, + h(x, ¥,))]. 
Thus 
PX Ver Vr 0) =F LOM(X rs Ve) + IM(Xe Ve) 
= MX Vr). 


The predictor-corrector method is consistent with the 
differential equation. 


5. Euler's method uses the recurrence relation 


Yui=¥+hY. 
For the differential equation 
y’ = —30y 


with »(0) = 1 and h = 0.1 we have 
Y41 = ¥,— 30 x O1Y, 
= -2Y,. 


With Yo = 1 we have the particular solution for this 
recurrence relation as 


¥ = (-2)". 


This is an increasing oscillating solution, (The sequence 
generated is 1, —2, 4, —8, 16, —32, 64...) 


The true solution is y = e~*°* which tends rapidly to zero. 
Something has clearly gone wrong with the method. 


6. Euler’s method is absolutely stable if —2 < ha <0. 
Hence for the differential equation 


y= —100y 
a = ~100 giving 
—2< -100h <0 


i. 100h < 2, (The inequality —100h <0 is always satisfied 
for positive h.) 


ie. h < 0.02, 


7. The differential equation 


y= —xyt 


1+ x? 
is linear with I(x) = —x and k(x) = u 
ineal (x) = —x =e 


The condition for absolute stability in Euler’s method is that 
—2<I(x)h <0. 

The minimum value for I(x) is — 10 since 0 < x < 10. 

Hence we require 
—2< —10h. 

ie =10h<2 


h<02, 


8, The theory tells us that y lies between 1000 and 2000 in 
the solution of the logistic equation 


y 11-2) with (0) = 2000, 


Hence the minimum value of 


dm y, 
rie 


occurs when y = 2000 giving $ = -30, 


We require 
dm 
—2enT <0 
ie 30h <2 
1 
easy 
<15 


Thus, provided h < 1/15, Euler's method is absolutely stable. 
9, This is a differential equation of the form 

y= m(x,y). 
AL x, we write 


MAY) = (Xp Sr): 


6) 
In this case m(x,y) = = 
so that 
y= 
m,(y) =——- 
y y 
Now 


The condition for absolute stability in Euler's method is 


dm, 

-2 <hr <0. 
x lies between 0 and 2 whilst y lies between 1 and 5, Hence 
the minimum value of 2 occurs when x, =2and y = 1 
giving 

dm, 


12 


Hence the method is certainly absolutely stable if 
=2< 12h, 
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1 
ie. h<=. 
ie h<e 


(In reality we could get a much less restrictive bound on h by 
obtaining the solution using the separation of variables 
method and then finding the minimum value of —6x/y*, but 
that would defeat the object of the exercise!) 


Solutions to the exercises in Section 4 
1. For the Taylor series method of order 2 we have 


he 


Yur=Y¥,+h¥e+>Yv. (1) 


Applied to the problem 


y) 

we have Y} = 4Y, and Y) = aY; =27Y, (differentiating the 
differential equation). Substituting into Equation (1) above 
gives 


2g? 
Yor= ¥, + hay, +2, 


ha? 
=(1+ha+——1Y,. 
( a + 2 
This is a recurrence relation of the form 
Y+1=aY, 
where 
Wo? 
=1+he+—. 
a=1+ha+ 2 


Note that this is the same equation for a that we had in the 
predictor-corrector method in tape frame 3. Hence the graph 
will be the same as in tape frame 4 and consequently the 
interval of absolute stability (|a| < 1) is (—2,0), 


In the logistic equation, with y(0) = 500 the solution 
increases rapidly to 1000 and 5 lies between 0 and —10 
where 


y 
= 10y{1 — ; 

m(y) (1 i dal 

As the interval of absolute stability for this method is (—2,0) 

the method is stable provided h > -2, ie. h< 02. 


2. Using the procedure for Simpson's method applied to 
linear differential equations with I(x) = 3, k(x) = sinx and 
h = 0.2 we obtain the recurrence relation 


Q2 .. F ; 
34 (sinx,-, + 4sinx, + sinx,.1) 


=Y¥,+ 15¥,-1+ Bling + 4sinx, + sinx,,;) 


Using the initial values ¥ = 0 at x =O and ¥, = 0.0255 at 
x, = 0.2 we have: 


Y, = 0.0255 + 15 x 0+ psino + 4sin0.2 + sin0,4) 


= 0.12417464; 
Ys = 0.12417464 + 1.5 x 0.0255 


a an 0.2 + 4sin0.4 + sin 0.6) 


= 0.35584007. 
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3. Using the procedure for linear equations with I(x) = 
k(x) = 0 and h = 0.1 we obtain the recurrence relation as 


y,, = [4x01 x10), | (3401 x 10 
mr \3—0.1 x 10 3—0.1 x 10, 


= 2Y,+2Y,-;. 
This is the required recurrence relation. 


¥,-1+0 


The auxiliary equation (cf. Unit 1, Section 2) is 
x? -2x-2=0. 
The roots of this quadratic are 


24+ /44+8 
2 


w=2x+2 0 or 


Am= 


=1+/3. 


7320508 and 
—0.73205081, 


The general solution for this recurrence relation is 
Yy=AKIM4+ Bx yh 

where A = 2.7320508 and 4 = —0.73205081, 

4. Again using the procedure for linear equations with 


I(x) = —10, k(x) = 0 and h = 0.1 we have 
, Bee 
mer \3= 01 x (=10)) 
3+ 0.1 x (10) 
+\3—01x (=10)) ot? 
=-Y¥+4Y,-1. 


This is a linear, constant-coefficient, homogeneous, second- 


order recurrence relation. Its auxiliary equation is 
wa=—-x¢h or x?+x-$=0. 


The roots of this equation are 


-livVl+2 
ie, d= — +2 « 0.3660254 and 


1,3660254, 


Mm 


2 
Hence the general solution is 
= A(0.3660254)" + B(—1.3660254)", 

If Y = 1, and Y, = 0.3679 the particular solution is obtained 
from solving for A and B in 

Y=A+B=1 

Y, = 0.3660254A — 1.3660254B = 0.3679. 
This gives A = 1.0010823 and B = —0,0010823. 
Hence the particular solution is 


¥, = 1,0010823 x (0,3660254)" 
—0.0010823 x (— 1.360254)", 


Compare this with the true solution 
Jy = (0.36787944)". 


The first term in Y, is almost exactly the same as the true 
solution while the spurious solution has a small coefficient 
(—0.0010823). However see what happens as n increases: 
(0.3660254)" tends to zero while (— 1.3660254)" is oscillating 
and increasing. Inevitably the spurious solution will 
eventually dominate the values obtained for Y,. Below are 
the values I obtained using the recurrence relation with 

Yo = Land Y, = 0.3679. 


r 0 1 2 3 4 5 


ae 0 0.1 02 0.3 0.4 0.5 


Y, 1 0.3679 | 0.1321 | 0.05185 | 0.0142 |0.011725 
; | 6 7 8 9 10 
| 0.6 0.7 08 09 1.0 


% {0.004625 10.0104875 |—0.0128}0.01804375 |—-0.02444375 


Note that the solution decreases initially, but then the 
spurious solution takes over so that after x = 0.5 the solution 
is oscillating and increasing. 


Solutions to the end of unit exercises in 
Section 6 


1. The recurrence relation for the trapezoidal method is 


Yor= Y,+ 500+ Yrea) 


h 
Y,+5(% + Yoo) + H2x, + 2x1) 


Hence 
h h 
Yall = ‘) = v(1 + ‘) + h(x, + Xy41): 


(Alternatively we could have used the procedure for linear 
equations with I(x) = 1 and k(x) = 2x.) 
With h = 0.1 we have 


LOSY, + O.1(x, + X41) 


2. The predictor-corrector method applied to the 
differential equation 


y'=Ty with y)=1 and h=0.1 
gives 

(i) Evaluate Y; = 7Y,. 

(ii) Predict Y#,, = ¥, +hY, = ¥, + ThY, =1.7¥,. 
(iii) Evaluate TY = 119Y,. 


(iv) Correct ¥,4; = Y,+ hay, + YF) 


=Y¥,+ toy, + 11.9Y,) = 1.945Y,. 


Hence Y,,; = 1.945Y,, 


The particular solution for this recurrence relation with 
Y% =1is 


¥, = (1.945)". 
3. For the recurrence relation 


h 
Yoi=¥, +54¥, + Visi) 


we have 


(rs Yes Yow ish) = 4(4¥, + Yes) 


Thus 
D(X VrsVrs0) = BAY + YF) 


= Vr = MX Vr) 


Hence the method is consistent with the differential equation. 


To compare this method with the trapezoidal method we 
compute the local truncation error. As it is an implicit 
method we can only compute the principal term, by 
assuming that all right-hand side values are correct; 


h 
ie. Veer = ye + 3(4y + rei) 
The Taylor series expansion for y.1 is 
re 
ort =e hye tye toe 
Hence 


3 
et 


The Taylor series expansion for y,«1 gives 


hk 
Yor = Yet hy + ove + 


h 
Voor =Ve thy, + 


Nees “Yeas = 7 Fp: 


3h? 
The principal term, — ay, involves h? while the principal 


term in the local truncation error for the trapezoidal method 
involves h?. We would thus expect better results using the 
trapezoidal method, 


4. For the problem 
yu 


we have 


h 
Yor = i+ g4a¥, + a¥,+1), 
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Collecting up terms gives 


he 4h: 
Yeos(1 —] = v(1 +), 


y.,, = (+ 4/5) 
hal) 
For absolute stability we require 


1 + 4ha/5) 
ha/S 


<i 


h 
If ha is positive i + “| > \ “ "| and the method is 


absolutely unstable, 


5 3 
which case the method is absolutely stable. 


0 
If ho is negative i + =| | T as ‘| provided ha > = in 


Hence the method is absolutely stable if — £ <ha <0. 


5, Using the procedure for linear equations with I(x) = 3 
and k(x) = 2x we obtain the recurrence relation for 
Simpson's method applied to the linear equation as 


rain (lye (My. 
+ LS + 4x, + X41): 
31 —h) 
Since x,-, + X,4, = 2x, this can be simplified to 


4h) [i+h 4h 
tan lg] (8 ghee 
To get started with this second-order recurrence relation, we 
require a value for Y,. This can be obtained, for example, by 


using a one-step method with very small step-size to compute 
an accurate approximation to Y,. 
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